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1  Introduction 


There  is  a  great  deal  of  current  interest  in  the  study  of  nonlinear  propagation  of  light 
in  periodic  structures  (see  e.g.,  the  recent  reviews  [1,  2,  3]  and  the  references  therein). 
The  combination  of  the  Kerr  nonlinearity  of  the  fiber  and  the  periodic  variation  of 
the  refractive  index  along  the  length  of  the  fiber  (fiber  grating)  produces  a  rich 
variety  of  phenomena  (bistability,  multistability,  switching  behavior,  etc.)  that, 
together  with  the  recent  developments  in  grating  fabrication  techniques,  gives  to 
these  periodic  structures  a  very  promising  potential  for  use  as  components  in  all- 
optical  communication  systems. 

The  localized  structures  that  emerge  from  the  balance  between  the  nonlinearity 
and  the  medium  periodicity,  called  gap  solitons,  exhibit  very  small  formation  lengths, 
enabling  nonlinear  pulse  compression  and  soliton  dynamics  to  be  observed  on  length 
scales  on  the  order  of  centimeters.  But  possibly  the  most  striking  feature  of  gap 
solitons  is  the  ability  to  propagate  (in  theory)  at  any  speed  between  zero  and  the 
speed  of  light  in  the  fiber  without  the  grating.  This  possibility  of  trapping  light  (gap 
soliton  propagation  at  zero  speed)  is  currently  a  topic  of  great  interest  because  of 
its  future  applications  in  the  production  of  all  optical  buffers  and  storing  devices. 

The  theoretical  research  on  light  propagation  in  fiber  gratings  is  based  on  the 
analysis  and  numerical  simulations  of  the  so-called  nonlinear  coupled  mode  equations 
(NLCME).  This  system  of  equations  accounts  for  the  effect  of  propagation,  nonlin¬ 
earity  and  the  medium  periodicity,  and  can  be  derived  from  the  Maxwell-Lorentz 
equations  (MLE)  for  electromagnetic  wave  propagation  in  a  dielectric  medium  as¬ 
suming  that  the  solution  is  a  superposition  of  two  counterpropagating  (backward 
and  forward)  fields  whose  envelopes  are  slowly  modulated  on  space  and  time.  In  the 
derivation  of  the  NLCME,  the  material  dispersion  terms  are  neglected  because  (clue 
to  the  slow  envelope  description)  they  are  small  as  compared  with  the  transport 
terms.  But  this  is  a  singular  problem  and  the  dispersion  terms  can  give  rise  to  small 
dispersive  scales  that  develop  in  the  propagation  time  scale.  This  phenomenon  is  not 
captured  in  the  NLCME  formulation  and  its  correct  description  requires  to  consider 
small  material  dispersion  terms. 

The  main  purpose  of  this  research  project  is  to  validate  the  widely  used  NLCME. 
That  is,  to  determine  whether  the  solutions  of  the  NLCME  constitute  good  approx- 
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imations  of  the  solutions  of  the  original  system  (i.e.,  the  material  dispersion  terms 
can  be  safely  ignored)  or  whether  the  dispersion  terms  are  essential  and  the  NLCME 
should  not  be  used.  As  a  secondary  goal,  this  research  project  also  aims  to  describe 
(applying  numerical  simulation  techniques  to  the  NLCME  with  dispersion  and  to 
the  MLE)  the  solutions  that  appear  when  the  NLCME  are  no  longer  applicable 
because  the  system  develops  small  dispersive  scales. 

This  report  is  organized  as  follows.  In  the  next  Section  we  introduce  our  physical 
model  for  light  propagation  on  a  fiber  grating:  the  one  dimensional  MLE,  and  study 
its  basic  linear  light  propagation  characteristics.  The  derivation  of  the  NLCME  with 
material  dispersion  terms  from  the  MLE  is  carried  out  in  Section  3.  The  simplest 
uniform  solutions  of  the  NLCME,  namely,  the  continuous  wave  solutions,  are  in¬ 
troduced  in  Section  4  and  its  linear  stability  properties,  including  the  effect  of  the 
material  dispersion,  are  analyzed  in  Section  5,  where  some  numerical  simulations 
of  the  NLCME  with  small  dispersion  are  also  performed  to  validate  the  theoretical 
predictions.  In  Section  6  we  describe  a  family  of  solitary  wave  pulse-like  solutions 
of  the  NLCME  known  as  Gap  Solitons,  and  in  Section  7  we  study  its  stability  prop¬ 
erties  theoretically  (defining  and  computing  the  Evans  function)  and  numerically 
(performing  some  numerical  simulations  of  the  NLCME).  Numerical  simulations  of 
the  MLE  are  used  in  Section  8  to  carry  out  a  definitive  check  of  the  stability  results 
and  of  the  validity  of  the  envelope  equations  with  dispersion.  And  finally,  the  main 
results  of  this  report  are  summarized  in  Section  9  and  the  details  of  the  numerical 
methods  used  to  integrate  the  NLCME  with  dispersion  and  the  MLE  are  given  in 
the  Appendixes  A  and  B. 
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2  Maxwell-Lorentz  equations 


We  describe  the  propagation  of  light  in  a  fiber  with  a  periodic  grating  and  a  cubic 
nonlinearity  using  the  one-dimensional  Maxwell’s  equations  [4,  5]  for  the  electro¬ 
magnetic  fields  evolution  together  with  an  anharmonic  Lorentz  oscillator  model  for 
the  polarization  (see  e.g.  [6,  7]  and  references  therein), 


dB  _  dE 

dt  dx  ’ 
dD  _  dB 

^ 0  dt  dx  ’ 

D  =  e0E  +  P, 

d2P 

^p2~dp  +  ^  _  2Ancos(2W-\?))P  -  7-P3  =  eo \E. 


(1) 

(2) 

(3) 

(4) 


In  the  system  above,  the  electric  field  E,  the  magnetic  field  B,  the  dielectric  dis¬ 
placement  D  and  the  polarization  P  are  scalar  fields  that  depend  on  the  spatial 
variable  x  and  on  time  t.  Hq  and  e0  denote,  respectively,  the  permeability  and  the 
permittivity  of  the  vacuum.  The  characteristic  frequency  Vtp  accounts  for  the  non 
instantaneous  polarization  response  of  the  media,  An  and  Xg  represent  the  strength 
and  the  period  of  the  grating,  that  is,  the  strength  and  the  period  of  the  spatial 
periodic  variation  of  the  refractive  index  of  the  fiber  (An  measures  the  size  of  the 
nonuniformities  of  the  refraction  index  relative  to  its  mean  value  no,  see  Fig.  1),  y 
is  the  linear  polarizability  of  the  medium  (jIq  =  1  +  y)  and  7  >  0  is  the  coefficient 
of  the  nonlinear  Kerr  effect. 


n{x) 


\ 9 


n0 


x 


Figure  1:  One  dimensional  fiber  with  periodic  variation  of  the  refractive  index. 

In  order  to  simplify  subsequent  calculations  it  is  convenient  to  make  the  system 
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(l)-(4)  non  dimensional  using  the  following  rescalings: 


B  =  vWM  )B,  D  =  (l/yfi)D,  E  =  (l/y/(^y))E,  P  =  (l/yfy)P, 

x  =  {\g/n)x,  t  =  (. \g/cn)t , 

here  c2  =  l/(eo fJ>o)  is  the  vacuum  speed  of  light.  After  dropping  tildes,  the  nondi- 
mensional  Maxwell-Lorentz  equations  (MLE  hereafter)  can  be  written  as 

dB  _  dE 

dt  dx  ’ 
d D  _  dB 

dt  dx  ’ 

D  =  E  +  P, 

d2P 

u ’  +  (1  -  2 An  cos(2a;))P  -  P3  =  (raj,  -  1  )E, 

y  dt- 

where  the  dimensionless  finite  time  polarization  response  frequency  is  now  a;2  = 
sW(  c2 7T2).  Notice  that  we  have  rescaled  the  spatial  variable  x  to  make  wavenum¬ 
ber  of  the  grating  equal  to  2  (see  eq.  (8)).  As  we  will  see  below,  the  reason  for 
this  choice  is  that  the  wave  trains  that  resonate  with  the  grating  and  develop  along 
the  fiber  will  then  have  wavenumber  1.  We  have  also  used  the  rescaling  to  absorb 
the  vacuum  properties  eo  and  /io  (the  vacuum  speed  of  light  is  now  equal  to  1  in 
the  rescaled  variables)  and  the  nonlinear  coefficient  7.  As  a  result,  the  MLE  above 
depend  only  on  three  dimensionless  parameters:  ca2,  An  and  n q. 

We  are  interested  in  the  description  of  the  solutions  of  the  MLE  in  the  physically 
relevant  regime  of  small  grating  strength  (i.e.,  near  uniform  refractive  index  fiber) 
and  small  nonlinear  Kerr  effect,  that  is, 

|An|  -C  1  and  \E\,  \B\,  \D\,  \P\  1. 

To  achieve  this  goal  we  need  to  understand  first  the  linear  light  propagation  charac¬ 
teristics  in  a  uniform  fiber.  The  next  step,  the  weakly  nonlinear  analysis  of  the  light 
propagation  in  a  nearly  uniform  fiber  that  is  performed  in  the  next  section,  is  just 
a  small  perturbation  of  this  basic  configuration  of  uniform  fiber  without  nonlinear 
effects. 

This  linear  analysis  is  carried  out  by  looking  for  uniform  wavetrain  solutions  with 


(5) 

(6) 

(7) 

(8) 
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spatial  wavenumber  k  and  frequency  uk, 

E(x,  t ) 

B(x,  t ) 

D(x,  t ) 
k  P(x,t) 

in  the  linearized  version  of  the  MLE  without  grating  An  =  0  (here  c.c.  stands  for 
the  complex  conjugate).  This  produces  the  following  algebraic  eigenvalue  problem 

^ kBk  kEk 

^kBk  kBk 
Dk  =  Ek  +  Pk 
(1  -  c jk2/u2p)Pk  =  \Ek 

that  has  nontrivial  solutions  if  and  only  if  ujk  satisfies  the  following  fourth  order 
polynomial  equation, 

ut  ~  ul(k2  +  LJpno)  +  ^pk2  —  0,  (10) 

which,  for  n( |  >  1  [4,  5,  8],  has  four  real  roots  of  the  form 

u)k  =  ±J(k2  +  WpW §)/2  ±  sj(k2  +  cn2ng)2/4  -  c (11) 
with  associated  eigenvectors, 

Ek  ul 

Bk  kujk 

= 

Dk  k2 

Pk  k2  —  u>l 

The  four  branches  of  the  dispersion  relation  (11)  are  plotted  in  Fig.  2.  The 
k  — >  —k  symmetry  of  uok  comes  from  the  invariance  of  the  MLE  without  grating 
(An  =  0)  under  spatial  reflections  x  — >  —x,  and  the  u>k  — >  —  u>k  symmetry  is  due 
to  the  fact  that  the  MLE  are  also  invariant  under  reflection  in  time  t  — »  —  t  (time 
reversal  symmetry).  For  every  wavenumber  k  there  are  four  possible  propagative 
wavetrains  of  the  form  given  in  eq.  (9),  and  the  corresponding  wa vet, rains  with 
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wavenumbers  —k  are  the  complex  conjugates. 


Figure  2:  Sketch  of  the  dispersion  relation  (11). 

It  is  interesting  to  notice  that  the  two  branches  of  u ok  that  remain  in  between  +k 
and  —k  correspond  to  wa vet, rains  with  polarization  in  phase  with  the  electric  held, 
i.e.,  k2  —  >  0  in  (12),  while  for  the  other  two  branches  the  fields  are  in  opposite 

phase.  There  are  also  two  very  different  behaviors  for  large  wavenumbers:  one  is 
dominated  by  the  finite  time  polarization  response  of  the  medium,  c ok  —■ ►  ±cnp  as 
k  — >  ±oo,  and  the  other,  cvk  — >  ±k  as  k  — >•  ±oo,  corresponds  to  propagation  like  in 
the  vacuum,  without  polarization  effects. 
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3  Envelope  equation  description 

In  this  section  we  derive,  from  the  MLE  (5)- (8),  the  envelope  equations  that  de¬ 
scribe  the  weakly  nonlinear  propagation  of  light  wavepackets  in  a  fiber  with  periodic 
variation  of  the  refractive  index. 

The  small  nonuniformities  of  the  refractive  index,  An  <C  1,  and  the  effect  of  the 
small  nonlinearity,  \B\,  \D\,  |P|,  \E\  <C  1,  can  be  accounted  for  by  allowing  the  linear 
wavetrains  obtained  in  the  previous  section  (9)  to  be  slowly  modulated  in  space  and 
time 

E(x,  t)  Ek 

B(x,t)  Bk 

'  >  = 

D(x,t)  Dk 

P{x ,  t)  Pk 

/  V  / 

Here  A  is  the  small  amplitude  of  the  wavetrain  whose  temporal  and  spatial  varia¬ 
tions  take  place  in  a  much  longer  scale  than  that  of  the  basic  wavetrain  (which  is 
assumed  to  be  of  order  one,  i.e.,  k  ~  1  and  ujk  ~  1).  The  expression  above  represents 
a  wavepacket  composed  of  modes  with  wavenumbers  in  a  narrow  band  around  k  and 
frequencies  near  ojk.  Using  standard  multiple  scales  techniques,  a  single  amplitude 
equation  (or  envelope  equation)  can  be  obtained  for  the  slow  evolution  of  the  ampli¬ 
tude  of  the  envelope  of  the  wavetrain  A(x,t).  This  kind  of  derivation  for  nonlinear 
optics  problems  can  be  found  in,  e.g.  [6,  7,  8],  and  in  the  review  [9]  for  more  general 
pattern  forming  phenomena  in  physical  systems  of  different  nature.  The  method 
used  to  derive  the  amplitude  equations  can  be  summarized  as  follows:  insert  the 
expansion  (13)  into  the  MLE  (5)-(8)  and  force  the  resonant  terms  to  cancel  at  every 
order  to  ensure  that  the  solution  remain  bounded  in  the  short  time  scale;  the  desired 
amplitude  equation  is  precisely  this  cancellation  condition. 

If  we  take  (13)  into  the  MLE  (5)-(8),  then  the  first  order  contribution  of  the 
grating  term  is  proportional  to 

2Ancos(2x)(Aeikx+iu>kt  +  c.c)  =  A  n{ei2x  +  e~i2x)(Aeikx+i  Ukt  +  c.c).  (14) 

The  expression  above  is  composed  of  wavepackets  with  frequencies  ±c ok  and  wavenum¬ 
bers  ±k  ±  2,  and  gives  rise  to  resonant  terms  (i.e.,  wavetrains  that  satisfy  the  dis¬ 
persion  relation  (11))  only  for  k  =  ±1  (see  Fig.  2).  Therefore,  as  anticipated  in  the 


A(x,t)  eikx+iUkt  +  C.C.  +  ....  (13) 
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previous  section,  only  those  wavepackets  (13)  with  wavenumber  and  frequency 


k  =  ±1 


Uk  =  w  =  \  (1  +  uln20)/2  ±J(  1  +  u2n20)2/A  -  uj2 


are  affected  by  the  grating.  The  effect  of  the  grating  manifests  itself  as  a  spatial 
resonance  that  couples  both  wavepackets:  the  wavepacket  with  k  —  1  produces  a 
resonant  term  with  wavenumber  k  =  —  1  and  the  same  frequency  and  vice  versa 
(see  eq.  (14)).  Wavepackets  with  wavenumber  fc  ^  ±1  do  not  feel  the  grating  in 
first  approximation,  they  propagate  like  in  a  uniform  fiber  and  the  weakly  nonlin¬ 
ear  evolution  of  its  amplitude  is  then  given  by  the  standard  nonlinear  Schrodinger 
equation,  see  e.g.  [4,  5,  7,  8]. 

In  the  remaining  part  of  this  section  we  derive  the  amplitude  equations  for  the 
weakly  nonlinear  evolution  of  two  resonantly  interacting  wavepackets  with  wavenum¬ 
bers  close  to  k  =  ±1. 

We  first  eliminate  the  dielectric  displacement  D  and  the  magnetic  field  B  from  the 
MLE  (5)-(8),  and  rewrite  them  as  a  system  of  two  second  order  partial  differential 
equations  for  the  polarization  P  and  the  electric  field  E, 


d2(E  +  P )  _  d2E 

dt2  dx 2  ’ 

d2P  2, 


=  — o>p(l  -  2Ancos(2x))P  +  u2p{n20  -  1  )E  +  ujpP3 


And  we  then  look  for  solutions  of  the  system  above  that  are  composed,  at  first  order, 
of  two  counterpropagating,  slowly  modulated  wavetrains  with  wavenumbers  k  —  ±1 
and  frequency  u  (15) 


E(x,  t ) 
P(x,t) 


=  V0(A+(x,t)eix+iu}t  +  A~(x,t)e~ix+iuJt)  +  c.c.  +  . . . ,  (18) 


where,  according  to  (12),  the  eigenvector  V0  is  given  by 


and  the  weakly  nonlinear  level  of  this  approach  requires  essentially  that 
•••  <  \A±X\  <  \A±\  <  1^1  <  1,  •  •  •  < \Af\  -C  1^1  <  1  and  An  -C  1,  (19) 


that  is,  small  amplitudes  that  depend  slowly  on  space  and  time  and  small  grating 
strength. 

The  appropriate  expansions  for  the  solution  of  eqs.  (16)-(17)  and  the  amplitude 
equations  in  powers  of  the  small  quantities  An,  A±,  At,  A±x,-  ■  .(which  will  be  treated 
here  as  independent  parameters)  are  thus  of  the  form 


E(x,  t ) 
P(x,t) 

At 

Ai 


|  =  Vo  (A+eix+iu)t  +  A~e~ix+iu}t)  +  c.c.  + 

+vt  At  +  A~  +  vt  Atx  +  v2  Axx  +  . . . , 

=  at  A+  +  at  At  +  at  Atx  +  ... , 

=  ttg  A  +  a1  Ax  +  a2  Axx  +  . . . . 


(20) 

(21) 

(22) 


Inserting  the  expansions  above  into  eqs.  (16)-(17),  a  linear  nonhomogeneous  system 
is  obtained  for  the  contribution  of  each  order.  For  the  resonant  terms,  i.e.,  those 
proportional  to  e±lx±lujt ^  a  condition  must  be  satisfied  to  ensure  that  there  are  not 
secular  terms  in  the  short  scale.  In  other  words,  the  linear  problems  corresponding 
to  the  resonant  terms  are  singular  and  hence  a  solvability  condition  must  be  satisfied 
by  the  nonhomogeneous  part;  these  solvability  conditions  give  the  coefficients  of  the 
amplitude  equations. 

Notice  that  only  the  resonant  terms  contribute  to  the  amplitude  equations  and 
only  the  amplitude  equation  for  A+  has  to  be  calculated  because  the  corresponding 
equation  for  A~  can  be  obtained  by  simply  applying  the  symmetry 


x  — >  — x  A+  < — »  A  , 


(23) 


which  comes  from  the  spatial  reflection  symmetry  of  the  original  problem  (5)- (8). 

The  linear  terms  in  the  amplitude  equations  can  be  easily  anticipated  because 
they  correspond  to  the  Taylor  expansion  of  the  dispersion  relation  (11)  at  k  =  1  (see 

e-g-  [9]), 


^k\k=i  -  U))A+  + 


duk 


dk 


At 


.1  d2u4 


k= 1 


2  dk 2 


Atr  + 


fc=l 
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The  first  coefficient  vanishes  (see  eq.  (15))  and  the  second  and  third  coefficients  are, 
respectively,  the  group  velocity  and  the  material  dispersion,  and,  after  making  use 
of  eq.  (10),  can  be  written  as 


dujk 


Vn  = 


id,  = 


~ul) 


dk  \k=i 

.1  d2ujk 

l~2  ~dkF 


UJA-UJ2p 


.  1  u>3(u2  —  1)  (u2  —  u>2)  (3u>2  +  u>4) 


k= 1 


-1-- 

2 


(^4  -  uiy 


(24) 

(25) 


The  first  order,  resonant  contributions  of  the  grating  to  the  expansion  of  the 
solution  (18)  and  to  the  amplitude  equation  (21)  are  of  the  form 


W AnA  elx+iu}t  and  wAnA  , 


where  the  two  component  vector  W  is  given  by  the  following  linear,  singular  non- 
homogeneous  problem 


u2-l 

a;2 

W  =-uj2v 

0 

0 

Vo  +  2iujw 

1 

1 

_  K  - 

u<2  ~  up 

p 

0 

1 

0 

1 

This  system  can  be  solved  only  if  the  right  hand  side  is  orthogonal  to  the  solution 
of  the  adjoint  problem 

w2  J  ' 

This  solvability  condition  gives  the  value  of  the  coefficient  of  the  amplitude  equation 


w  = 


.  Ld(  1  —  U2) 

*2K-w2) 


(26) 


The  first  order  contributions  of  the  nonlinear  term, 


f/iA+IA+lV3^,  U2A+\A~\2eix+iuJt  and  UlA+\A+\2,  u2A+\A~\2, 


are  computed  similarly.  The  following  linear  problems  are  obtained  for  the  vectors 
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U\  and  U2 


u2  -  1  u2 

U\  —  — 3cu2 

0 

+  2itoui 

1  1 

_  (™o  -  1  Vp  _ 

1  P 

(1  —  cu2)3 

O 

1 _ 

cu2  —  1  cu2 

_  (no  -  1  Vp  -  u2v 

U2  =  -6c u2 

0 

(1-w2)3 

+  2iuu2 

1  1 

0  1 

and,  after  applying  the  solvability  condition,  the  resulting  amplitude  equation  coef¬ 
ficients  are  given  by 


u  1 


,3cu(l  —  uj2)3 
1  2(o;4  -  u® 


OL 


and  112 


,3cu(l  —  cu2)3  2 


(27) 


Collecting  the  coefficients  above  (24)-(27)  and  applying  the  spatial  reflection  sym¬ 
metry  (23)  the  amplitude  equations  can  be  finally  written  as 


A+  —  VgA+  +  idA+x  +  10  A nA  +  ^4^ (xti | ^4”*” | 2  T  U2\A  |2)  +  . . . ,  (28) 

At  =  — VgAx  +  idAxx  +  w/S.nA^  -\-  A  (ui|4  |2  +  1/2 1  A3- 12)  +  . . . .  (29) 


Notice  that  the  grating  and  nonlinear  coefficients  (26), (27)  are  purely  imaginary. 
This  result  could  have  been  anticipated  without  any  calculations  because  the  am¬ 
plitude  equations  have  to  be  invariant  under  the  transformation 

t  ->  -t  A±  -»•  (30) 


that  comes  from  the  time  reversal  (t  — >  —t)  invariance  of  the  MLE  (16), (17).  Also, 
the  fact  that  the  only  nonlinearity  in  the  MLE  is  cubic  forces  the  nonlinear  coeffi¬ 
cients  (27)  to  verify 

U2  =  2u\. 

We  will  consider  the  simplest  possible  geometrical  configuration:  propagation  of 
light  in  a  fiber  ring  with  length  L  >  1.  The  spatial  periodicity  condition  for  the 
electric  held  and  the  polarization  implies  that  the  boundary  conditions  for  A+  and 
A~  are  (see  eq.  (18)) 

A+(x  +  L)el6  =  A+(x,  t),  A~(x  +  L)e~l6  =  A~(x,  t).  (31) 
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Here  6  =  L  (mod27r)  measures  the  mismatch  between  the  natural  wavelength  of  the 
wavepackets  (=27r)  and  the  period  of  the  domain,  but  we  will  restrict  our  analysis  to 
the  particular  case  6  =  0  (i.e.,  ring  length  equals  to  an  integer  multiple  of  the  basic 
period  of  the  wavepackets)  to  reduce  the  number  of  parameters  of  the  problem. 

There  are  two  possible  choices  for  u  (o^)  depending  on  the  sign  selected  in  (15), 
see  Fig.  3.  As  mentioned  in  the  previous  section,  these  two  possibilities  correspond  to 
a  configuration  with  the  electric  field  and  the  polarization  in  phase  (for  u  =  u~)  and 
in  opposite  phase  (for  u  =  cu+).  The  group  velocity  (24)  is  positive  in  both  cases  (it 
is  the  slope  of  the  curve  04.  at  k  =  1  in  Fig.  3)  but  the  sign  of  the  material  dispersion 
coefficient  d  (25),  which  is  related  to  the  curvature  of  the  curve  04  in  Fig.  3,  changes. 
The  nonlinear  and  grating  terms  have  imaginary  parts  that  are  always  negative;  see 
eqs.  (26)  and  (27)  and  Fig.  3,  and  recall  that,  using  the  dispersion  relation  (10)  for 
k  =  1,  the  denominator  can  be  written  as  u>A  —  uj2  =  u2[(cj2  —  1)  +  (t J2  —  uj2n^)\. 


Figure  3:  Detail  of  the  dispersion  relation  (11)  with  the  two  frequencies  uA  for  k  =  1. 

In  order  to  make  the  nonlinear  and  grating  coefficients  positive,  we  will  work  with 
the  complex  conjugates  of  the  amplitudes  and,  to  absorb  some  of  the  parameters  of 
the  problem,  we  will  also  perform  the  following  rescalings 

x  =  Lx,  t  =  ( L/vg)t ,  A±  =  \Jvg/(L\u2\)A±,  (32) 
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that,  after  dropping  tildes,  yield  the  scaled  equations 


(33) 

(34) 

(35) 


At  =  AZ  +>£41+  +  mA~  +  i4>|.4+|2  +  |  A'  |2), 
A r  =  -A;  +  isA~x  +  inA+  +  (cr|^- 12  +  |Vf), 

A±(x  +  1,  t)  —  A±(x,  t), 


where  £  =  — d/(Lvg )  C  1  is  positive  (negative)  for  u>  =  oj+{oj  =  u  ),  the  scaled 
grating  strength  k,  =  AnL\w\/vg  ~  1  is  always  positive,  and  the  nonlinear  coefficient 

*=§• 

If  we  neglect  the  material  dispersion  terms  (i.e. ,  set  £  =  0)  in  the  system  above 
then  the  so-called  nonlinear  coupled  mode  equations  (NLCME)  are  obtained, 


A+  -A+  =  mA~  +iA+(cr\A+\2  +  \A~\ 2),  (36) 

A~  +  A~  =  inA+  +  iA~(a\A~\2  +  |A+|2),  (37) 

A±(x  +  1,  t)  —  A±(x,  t),  (38) 


Equations  (36)-(37)  represent  a  balance  of  the  transport  at  the  group  velocity,  the 
effect  of  the  grating  and  the  Kerr  nonlinearity.  This  hyperbolic  system  and  its  more 
remarkable  solutions  have  been  widely  studied  theoretically  and  numerically,  see  e.g. 
[10,  11,  2,  6,  12,  13,  15,  16,  17,  18,  19,  20,  14].  Also,  some  interesting  generalizations 
of  the  NLCME  have  been  recently  analyzed,  including  dual  core  fibers  and  quadratic 
nonlinearities  [3],  stimulated  Brillouin  scattering  [21]  and  nonlinear  effects  in  the 
grating  strength  [22,  23].  And  similar  hyperbolic  systems  have  been  used  for  the 
description  of  Bose-Einstein  condensates  [24]  and  spontaneous  Hopf  bifurcations  in 
dissipative  pattern  forming  systems  [25,  26,  27]. 

The  purpose  of  this  research  project  is  to  study  the  effect  of  the  addition  of  small 
dispersive  terms  to  the  solutions  of  the  NLCME,  that  is,  to  analyze  the  system  (33)- 
(35)  in  the  limit  £  — >  0.  There  are  some  recent  papers  [28,  29,  30]  that  study  the 
NLCME  with  dispersive  terms  but  they  do  not  cover  the  physically  relevant  regime: 
small  dispersion  coefficient. 

The  small  dispersive  terms  of  the  system  (33)-(35)  constitute  a  singular  pertur¬ 
bation  of  the  NLCME  and  introduce  a  new  small  dispersive  scale  Id  rs_/  i- 

The  solutions  of  the  NLCME  are  approximated  solutions  (up  to  0(e)  errors)  of  the 
system  with  small  dispersion  (33)- (35).  Our  goal  is  to  determine  which  solutions  of 
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the  NLCME  remain  stable  in  the  system  with  small  dispersion  and  which  ones  are 
unstable  against  perturbations  containing  small  dispersive  scales.  For  the  solutions 
that  develop  dispersive  scales  the  dynamics  predicted  by  the  NLCME  is  not  correct 
and  the  complete  system  with  dispersion  must  be  integrated.  Similar  analyses  can 
be  found  in  [27]  in  the  context  of  the  oscillatory  instability  in  dissipative  systems, 
and  in  [31]  for  the  onset  of  the  Faraday  instability  in  a  nearly  conservative  system. 

The  next  sections  are  dedicated  to  the  study  the  effect  of  the  small  dispersion 
in  a  well  known  family  of  solutions  of  the  NLCME,  namely,  the  continuous  wave 
solutions. 
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4  Continuous  wave  solutions 

The  simplest  solutions  of  the  NLCME  (36)- (38)  are  the  so-called  continuous  wave 
fields  (CW)  that  correspond  to  amplitudes  with  steady,  uniform  modulus  (see  [17] 
and  [20]).  These  solutions  are  of  the  form 

A+  =  A+W  =  p  cos  9  eiat+imx,  (39) 

A~  =  A~w  =  p  sin  9  eiat+imx ,  (40) 


where  p  >  0  and  9  e  [— f ,  f],  and  the  frequency  and  wavenumber  of  the  amplitudes 
are  given  by 


K 


a  = 


cr  +  1 


sin  29 


-P‘, 


m  =  (  — 


K 


a 


sin  2  9 


~-p2)  cos  2 9. 


(41) 

(42) 


The  corresponding  physical  fields  (see  eq.  (18))  are  the  superposition  of  two  uniform 
counterpropagating  wavetrains,  and  m  and  a  represent,  respectively,  small  correc¬ 
tions  to  the  basic  wavetrains  wavenumber  and  frequency.  To  be  more  precise,  the 
resulting  wavenumbers  and  frequencies  of  the  wavetrains  in  eq.  (18)  are  1  ±  ^  and 
uj+  fvg,  with  L  1.  The  parameter  p  represents  the  total  power  propagating  along 
the  fiber, 

P2  =  \Atw\ 2  +  IA:J2> 


and  9  measures  the  ratio  between  the  two  wavetrains.  For  9  =  ±|  both  amplitudes 
have  the  same  size  and  the  physical  fields  take  the  form  of  a  standing  wave,  while 
for  9  — ■>  0  and  9  — >  one  of  the  amplitudes  vanishes  and  the  pattern  inside  the 

fiber  is  a  pure  travelling  wave. 

In  the  linear  limit  of  small  light  intensity  (p  =  0)  the  relation  between  the  CW 
wavenumber,  m,  and  frequency,  a,  is  given  by 

2  2/i  \ 

">  =  a  1 - j  , 

az 

which  is  plotted  in  Fig.  4.  It  is  clear  form  this  figure  that  the  main  effect  of 
the  grating  in  the  linear  light  propagation  characteristics  (discussed  at  the  end  of 
Section  2,  see  Fig.  2)  is  to  open  a  frequency  gap  of  size  2k,  centered  at  the  resonant 
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point  ( m,a )  =  (0,0),  i.e.,  (k,Uk)  =  (1,cp)  in  terms  of  the  original  variables  (9), 
for  which  there  is  no  possible  propagation.  This  opening  of  a  “photonic  band  gap” 
of  excluded  frequencies  is  a  typical  resonance  effect  in  wave  propagation  problems 
in  linear  periodic  media.  For  the  linear  modes  with  wavenumbers  near  rn  —  0 
(9  close  to  in  eqs.  (41)  and  (42),  see  Fig.  4)  the  reflection  produced  by  the 
grating  is  maximum  and  when  we  move  away  from  the  resonance  point  (0^0  and 
6  — >  ±| )  the  effect  of  the  grating  is  gradually  reduced  and  the  linear  propagation 
characteristics  of  bare  fiber  (a  =  ±m,  in  our  scaled  variables)  are  recovered. 


Figure  4:  Linear  light  propagation  characteristics  in  a  fiber  grating  showing  the  frequency  gap. 


As  the  power  is  increased,  p  >  0,  the  effect  of  the  nonlinearity  of  the  fiber  comes 
into  play  and  the  relation  between  the  frequency  and  the  wavenumber  of  the  CW 
becomes  more  complicated, 


m2  =  (a  —  up  )"(1 


K 


[a 


'  1 +<7 
-  2 


)p 


2\2  ' 


The  equation  above  is  a  fourth  order  polynomial  equation  in  a  and  it  does  not  have 
real  solutions  if 

\a  —  ( — 2 — )p  I  - 

hence  the  frequency  gap  remains  of  size  2k  and  it  is  just  shifted  upwards  by  the 
nonlinear  effects.  The  structure  of  the  CW  is  essentially  similar  to  that  of  the  linear 
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case,  i.e.,  there  are  two  CW  with  for  every  wavenumber  m,  if  the  power  is  below  the 
critical  value 


Pc  = 


5 


see  Fig.  5a,  while  for  p  >  pc  the  nonlinear  distortion  becomes  more  important  and 
a  region  of  higher  multiplicity  of  CW  solutions  (up  to  4)  develops  near  rn  —  0  (see 
Fig.  5b). 


Figure  5:  a  —  m  plot  of  the  CW  with  constant  power  p  for  a)  p  <  pc  and  b)  p  >  pc  (a  —  |). 


We  conclude  the  description  of  the  CW  introducing  the  representation  of  this 
family  of  solutions  in  the  ( 6,p 2)  plane,  shown  in  Fig.  6.  This  is  a  convenient  way 
to  represent  the  CW  because  each  point  on  the  ( 9 ,  p2)  plane  corresponds  to  a  single 
CW  and  vice  versa,  and  there  are  not  multiple  solutions.  In  Fig.  6,  white  (shaded) 
regions  correspond  to  positive  (negative)  wavenumber,  and  the  thick  lines  are  the 
rn  —  0  curves,  that  are  given  by,  see  (42), 


p2  =  pI  = 


2k 


(1  —  a)  sin(2d) 


and  6  =  ±2 


(43) 


The  CW  inside  the  region  p2  >  p\  in  Fig.  6  correspond  to  those  that  lie  along  the 
loop  that  appears  for  p  >  pc  in  the  a  —  m  plots  of  the  CW,  see  Fig.  5b.  Constant 
m  (m  7^  0)  lines  are  plotted  in  Fig.  6  as  thin  lines  with  arrows  indicating  the  m 
growing  direction  (recall  that  m  — >  ±oo  for  6  — >  0  and  6  — >  ±f ,  see  eq.  (42)).  Note 
that  the  only  allowed  values  of  the  wavenumber  m  are  those  compatible  with  the 
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periodicity  boundary  conditions  (38),  i.e.,  those  of  the  form 


m  =  27m,  with  n  e  Z. 


2  4  4  2 


Figure  6:  The  ( 6,p 2)  representation  of  the  CW  (a  =  |). 

It  is  interesting  to  notice  also  that  this  ( 9 ,  p2)  representation  is  somewhat  redun¬ 
dant  because  the  CW  corresponding  to  9  and  that  corresponding  to  ±|  —  9  are  the 
same  after  a  spatial  reflection  (23).  This  results  in  the  symmetry  around  9  —  ± 
that  can  be  appreciated  in  Fig.  6. 
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5  Stability  of  the  CW 

The  linear  stability  characteristics  of  the  family  of  continuous  wave  (CW)  solutions 
(39)-(40)  are  obtained  in  this  section.  This  is  performed  in  two  steps:  we  will  first 
study  the  stability  of  the  CW  against  perturbations  with  wavenumber  of  order  unity 
(this  completes  the  analysis  presented  in  [17]),  and  we  will  then  use  the  NLCME 
with  small  dispersion  terms  (33)- (35)  to  compute  the  stability  of  the  perturbations 
that  contain  small  dispersive  scales. 

The  evolution  of  the  infinitesimal  perturbations  of  the  CW,  defined  by 

A+  =  A+W(l  +  a+),  AT  =  A~w(l  +  a-),  with  |a±|  <C  1, 

is  given  by  the  linearized  version  of  equations  (33)-(35) 

a~t  —  ax  =  in{a~ —  a+ )  ta nd+icrp2  cos2  9(a++  a,+)+jp2  sin2  9(a~+  a~)+i£axx, 
a~[  +  a~=  in{a+—  a~ )/ ta nO+icrp2  sin2  9(a~+  a~)+ip2  cos2  9(a++  a+)+iea~x, 
a±{x  +  1  ,t)  —  a±(^,  t). 


This  linear  system  is  solved  via  the  Fourier  expansion 

OO 

(a+,a-)=  ]T  «(t).Oi(t))ei2'fa, 

k=—oo 

where  k  represents  the  spatial  wavenumber  of  the  perturbation  and  the  coefficients 
satisfy  the  following  system  of  ordinary  differential  equations 


i(2irk)a k  +in(ak  —  a k)  ta nO+icrp2  cos2  9 (a k  +  a+k) 
+ip2  sin2  9(ak  +  a~k)  -  i£(27rk)2ak , 

— i(27rk)ak  +in(ak  —  ak)/  ta  n9+iap2  sin2  9{a] 7+  a~k) 
+ip 2  cos2  9(ak  +  atk)  —  i£(2nk)2ak  . 


(44) 

(45) 


The  small  parameter  £  in  the  system  above  allows  us  to  distinguish  between  two 
essentially  different  types  of  perturbations  that  are  studied  separately:  perturbations 
with  (k  ~  (1  / \f\£~\)  3>  1)  and  without  (k  ~  1)  small  dispersive  scales. 
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5.1  Perturbations  with  k  ~  1  wavenumber 


In  this  case  the  dispersive  terms  in  eqs.  (44)- (45)  are  small  and  can  be  neglected. 
The  resulting  equations  together  with  those  corresponding  to  a±k  form  a  quartet 
that  is  uncoupled  from  the  equations  for  the  other  wavenumbers.  If  we  look  for 
solutions  of  this  quartet  proportional  to  e  ,  the  following  fourth  order  polynomial 
equation  for  12  is  obtained 

(122  +  (27t k)2)2  +  2k2(122  +  ( 2nk )2)  +  4 np2  t&n^2  ((27r/c)2(l  +  a)  +  122(1  -  a)) 

1  +  tan  9 

K2 

+k2  tan2  0(12  +  i(27tk))2  H - 5— -(12  —  i(2irk))2  =  0, 

tan  0 

that  can  be  turn  into  a  real  coefficient  polynomial  equation  by  means  of  the  change 
of  variable  12  =  iuj 

(t u2  —  (2nk)2)2  —  2k2(u2  —  (2irk)2)  +  4  up2 - ((27tA;)2(1  +  a)  —  cu2(l  —  a)) 

1  +  tan  0 

—k2  tan2  9(lu  +  (2nk))2 - ^-(w  —  (27tA;))2  =  0.  (46) 

tan  0 

For  the  particular  case  of  uniform  perturbations,  k  =  0,  the  solutions  of  (46)  can 
be  calculated  explicitly 


c 0  = 


± 


u  =  0  (double) 


4k2 


sin2  (29) 


2/qo2(l 


and 


—  a)  sin(20), 


and  thus  we  have  instability,  i.e.,  12  with  positive  real  part,  when 


2  .  2 

P  >  Po 


2k 

(1  -a)  sin3  (29)' 


(47) 


This  instability  region  lies  inside  the  higher  multiplicity  region  p2  >  p2,  see  eq. 

(43),  and  is  represented  in  Fig.  7.  The  condition  above  states  that  the  CW  that 

correspond  to  the  upper  part  of  the  loop  represented  in  Fig.  5b  are  unstable.  This 

can  be  easily  seen  from  expression  (42):  the  vertical  slope  points  in  the  a  versus  m 

din 

curve  in  Fig.  5b  are  given  by  the  condition  -7^—  =  0  that  is  equivalent  to  p2  —  p2y 
In  the  limit  of  large  wavenumbers,  k  ±00,  the  solutions  of  eq.  (46)  can  be 
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expanded  as 

lu  =  lu0  +  lui  +  ...  i  with  |cu0|  S>  |cui|  3> - 

The  first  order  term  cu0  ~  k  3>  1  satisfies  the  equation 


K2  -  (2 nk)2)2  =  0, 


whose  solutions  are  of  the  form  cuo=±27tA;  and  correspond  to  purely  imaginary  values 
of  12.  In  order  to  determine  the  stability  of  the  CW  we  have  to  compute  the  next 
order  correction  u\  ~  1  that  is  given  by 


9  9  _i_o  „  9  tan  9 

cut  =  k  tan  9  —  2a up  - 77—. 

1  r  1  +  tan2  9 

The  CW  are  then  unstable  if  the  imaginary  part  of  cui  is  nonzero,  and  this  occurs 
inside  the  region  defined  by  the  conditions 


2  ^  2 

P  >  Poo 


k tan±2 9 
a  sin (29)  ’ 


(48) 


that  is  represented  in  Fig.  7.  Notice  that,  for  the  CW  in  this  region,  all  perturbations 
with  wavenumber  above  a  certain  threshold  are  unstable.  This  instability  affects  only 
the  upper  branch  of  CW  in  Fig.  5a  ( 9  >  0)  and  if  the  power  on  the  fiber  is  such 
that  p 2  >  k/a  then  all  CW  along  this  branch  are  unstable. 

The  stability  analysis  of  the  CW  is  completed  with  the  determination  of  the 
unstable  perturbations  with  finite  nonzero  wavenumber.  These  perturbations  cor¬ 
respond  to  complex  conjugate  pairs  of  roots  (with  nonzero  imaginary  part)  in  eq. 
(46), 

cu  =  a  ±  ib  {b  ^  0)  12  =  ±6  +  ia  unstable. 


This  task  can  be  somewhat  simplified  by  taking  into  account  the  invariance  of  eq. 
(46)  under  the  transformations 


(cu,  kP)  — >  ( — cu,  — kP)  and  (cu,  tanf?)  — 4  (— u,  1/ tand) 


that  allows  us  to  restrict  the  search  to  the  parameter  range 

_  1  _  ,.7171", 

k  >  0  and  9  e - ,  —  . 

4  4 
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Figure  7:  Stability  properties  of  the  CW  for  er  =  ^  (shading  indicates  instability).  Thin  solid  lines 

correspond  to  k\  =  0.25  +  0.25n,  with  n  =  0, 1,2, _  For  a  given  value  of  k,  the  CW  above  the 

line  k\  =  2tt/k  are  unstable. 

Furthermore,  if  we  use  the  scaled  variables 

Cu  =  u/k,  k=(2nk)/hi  and  p2  =  p2/n,  (49) 

then  the  strength  of  the  grating,  k  >  0,  can  be  absorbed  and  eq.  (46)  can  be 
rewritten  as 

(<52  -  k2)2  ~  2(d2  -  k?)  +  4p2_J|4_(^(1  +  a)  -  i>2(l  -  a)) 

—  tan2  8(ll>  +  k)2 - (o’  —  k)2  —  0. 

tan  0 
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The  onset  of  complex  roots  takes  place  when  both  the  equation  above  and  its 
derivative  with  respect  to  Co  vanish.  From  these  two  equations  the  values  of  9  and  p2 
are  numerically  computed  for  any  given  values  of  Co  and  k,  and  the  following  results 
are  obtained: 

•  If  6  >  0  then  there  are  no  more  instabilities  apart  from  the  one  already  obtained 
in  the  limit  k  — >  ±oo  (48),  see  Fig.  7. 

•  For  negative  6  and  for  each  value  of  p2  =  p2  j k  there  is  a  critical  scaled  wavenum¬ 
ber  k\  such  that  the  perturbations  with  k  <  k\  are  unstable  and  those  with 
k  >  k\  are  stable  (the  curves  of  constant  k\  value  are  plotted  as  thin  solid  lines 
in  Fig.  7).  A  CW  will  be  then  unstable  if  the  perturbation  mode  with  lowest 
wavenumber  ( k  =  2tc)  is  unstable,  that  is,  if 

97 r 

—  <  ku  (50) 

K 

see  (49).  Note  that  the  previous  instability  boundaries  (47)  and  (48)  depend 
only  on  the  combination  of  parameters  p2  /  n,  while  the  new  instability  defined 
by  the  condition  above  depends  also  on  k;  for  a  given  value  of  k,  all  the  CW 
above  the  line  k\  =  / n  are  unstable,  see  Fig.  7. 

•  The  CW  near  the  p2  =  p$  line  present  a  small  interval  of  stable  scaled  wavenum¬ 
bers  k  G  [0,  k2\,  but  this  does  not  modify  the  instability  region  defined  above 
because,  as  it  has  been  checked  numerically,  k2  remains  always  below  k\/2. 

Therefore,  when  we  enter  the  region  defined  by  condition  (50)  the  CW  become 
unstable  against  perturbations  with  wavenumber  k  =  ±1,  that  is,  with  one  single 
wavelength  along  the  fiber.  The  lines  of  constant  value  of  k\  mark  another  instability 
limit  for  the  CW,  these  lines  move  towards  higher  values  p2  /  n  as  k\  increases  (see 
Fig. 7)  and  they  enter  the  higher  multiplicity  region  p2  >  p2x,  see  eq.  (43),  at  the 
critical  value  kic  =  2A/(T+_o:)7(T^o:). 

We  can  summarize  the  results  on  this  section  as  follows.  The  linear  stability  of  the 
CW  (without  considering  dispersion  effects)  is  determined  by  the  three  conditions 
given  by  eqs.  (47), (48)  and  (50),  and,  depending  on  the  strength  of  the  grating  k, 
we  can  have  two  different  configurations  that  are  sketched  in  Figs.  8  and  9.  The  CW 
with  positive  9,  i.e.,  those  on  the  upper  branch  of  the  a  —  m  representation  in  Fig.  5, 
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are  always  destabilized  via  the  high  wavenumber  instability  (48).  While  for  negative 
6  the  situation  is  more  complicated:  if  k  >  nc  =  2n/kic  =  7 Ty/ (1  —  <r)/(l  +  cr)  then 
the  hnite  wavenumber  instability  (50)  dominates  (Fig.  8),  but  for  k  <  nc  the  zero 
wavenumber  instability  (47)  is  also  present  and  destabilizes  the  CW  with  6  near  —j, 
see  Fig.  9.  Notice  that  for  large  values  of  the  power  in  the  grating  p 2  all  the  CW 
are  unstable  and  that,  as  the  power  is  increased  the  first  CW  to  become  unstable 
are  those  that  are  close  to  pure  TW,  that  is,  those  with  6  near  0  or  ±7t/2. 


Figure  8:  Stability  properties  of  the  CW  for  cr  =  ^  and  n  >  kc  (shading  indicates  instability) 
together  with  the  corresponding  a  —  m  plots  (see  Fig.  5)  for  several  representative  values  of  p2 
(solid  and  dashed  lines  indicate  stable  and  unstable  CW,  resp.). 


5.2  Perturbations  with  small  dispersive  scales 

Figs.  8  and  9  summarize  the  stability  properties  of  the  CW  in  the  context  of  the 
NLCME  (36)- (38),  that  is,  without  small  dispersion  terms  (e  =  0  in  eqs.  (33)- 
(35)).  In  the  limit  £  — >  0,  the  effect  of  the  dispersion  terms  on  the  stability  of 
the  CW  amounts  to  just  a  small  correction  of  the  results  obtained  in  the  previous 
section  (e  =  0)  for  perturbations  with  wavenumber  k  ~  1.  This  is  not  true  for  the 
perturbations  with  wavenumber  on  the  dispersive  scale,  i.e.  \k\  i/7RT»K  see 


24 


l 

a 

X.  m 

1 

v 

a 

)  r 

m 

Figure  9:  Stability  properties  of  the  CW  for  a  =  \  and  n  <  kc  (shading  indicates  instability) 
together  with  the  corresponding  a  —  to  plots  (see  Fig.  5)  for  several  representative  values  of  p 2 
(solid  and  dashed  lines  indicate  stable  and  unstable  CW,  resp.). 


eqs.  (44)-(45)).  For  these  perturbations  the  dispersive  terms  must  be  retained  and 
their  stability  characteristics  can  be  calculated  using  multiple  scales  techniques. 

To  this  end  we  define  the  scaled  wavenumber  K  =  (2nk)^/\£\  ~  1  and  expand 
the  solution  of  eqs.  (44)-(45)  in  the  form 


aK  ~  aK0 

where  T  = 


(f,  T )  +  \f\e\  T)  H - , 


K 


t/ y/\e\.  At  first  order  we  obtain 


^k-iKa 
dT 

da 


KO 


K  0 


dT 


+  iKa 


KO 


—  aKo(ti  T)  +  y/\e\  aK1(t,  T)  H - 

=  0, 

=  0, 


5 


whose  general  solution  is  given  by 


1K01  a  K  0 


3)  ~  (^A'o(^)e 


iKT  A- 


)(t)e~ 


A-  ( „—iKT\ 
AK  O'- 


Thus,  in  the  fast  time  scale  T,  the  group  velocity  term  dominates  and  these  short 
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wave,  dispersive  perturbations  simply  travel  in  opposite  directions.  At  the  next 
order  we  obtain 


pp  -  iKa% i  =  [ — -  i(K  tan#  +  j^K2)A+0  +  iap 2  cos2  #(A+0  +  A±A0)]elKT 

+[«  tan  0Ako  +  ip2  sin2  #(Aao  +  AZK0)\e~lKT , 

pp  +  iK a~Kl  =  [-^p  -  *(k/  tan  0  +  pAT2)AA0  +  iap2  sin2  #(Aao  +  AlA0)]e"^T 

+[ik/  tan  #Aao  +  ip 2  cos2  #(Aao  +  A+A0)]eiKT. 

The  higher  order  correction  (aA1,aA1)  will  then  remain  bounded  in  the  fast  time 
scale  T  if  the  following  equations  are  satisfied 

pjp  =  ~i{K tan6>  +  -^K2)A^0  +  iap 2  cos2  9(A~^0  +  A+ao), 

(JA-  F  _ 

— p  =  — *(/c/  tan  6*  +  —  K' 2)Aj<0  +  iap2  sin2  #(Aao  +  AzK0). 

These  equations  together  with  the  corresponding  ones  for  (AlA0,  AZko)  form  a  linear 
system  with  constant  coefficients  that  give  the  evolution  the  dispersive  perturbations 
with  wavenumber  ±K  in  the  slow  time  scale  t.  The  left  (+)  and  right  (“)  propa¬ 
gating  perturbations  (associated  with  the  amplitudes  A+and  A~,  respectively)  are 
uncoupled  and  their  stability  properties  are  given  by  the  exponents 

=  ±A/ptan#  +  —  K 2){2ap2  cos2  6  —  (ft  tan  9+  —  A'2)),  (51) 

V  p|  \£\ 

=  ±P(k/  tan ~9  +  —  K2)(2ap2  sin2  9  —  (k/  tan  9  +  —A'2)),  (52) 

and  hence  are  unstable  if  the  following  conditions  are  satisfied 

0  <  k  tan  9  +  pA'2  <  2a p2  cos2  9, 

|e| 

0  <  k/  tan  9  +  pA'2  <  2a p2  sin2  9 , 

£ 
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that  can  be  expressed  as 


2  \  2 
P  >P+ 


2  \  2 
p  >p- 


tan  9  ,  „  e 

- — —  (Kt&nd  +  —  K1) 

a  sin (20) v  |e|  J 


£  K2 

with  tan  9  >  —  — ,  and 

£  K 


tan  1 9  ,  „  £  o i  ,  „  £  K2 

- — —  (k  tan  19+1—rK2)  with  tan  1  6  >  — -  — 

a  sin (20) v  |e|  ;  “  |e|  /c 


(53) 

(54) 


Notice  that  if  we  set  Jl 2  =  0  then  the  instability  condition  (48)  is  recovered,  in  other 

words,  the  result  for  the  limit  k  — »  ±oo  in  the  k  ~  1  regime  matches  with  the  result 

for  K  — >  0  in  the  k  ~  1/  ^/jej  S>  1  regime. 

The  above  instability  criterion  define  a  family  of  regions  in  the  6  —  p2  plane,  which 

depends  on  the  parameter  K 2.  The  expression  (53)  is  identical  to  (54)  if  we  change 

tan  9  by  1/  tan  9,  therefore  it  is  enough  to  obtain  the  instability  region  defined  by 

the  p2+  condition;  the  corresponding  one  for  p\  is  simply  the  symmetric  around  the 
7 r 

vertical  axes  9  —  ±— . 

4 

Fig.  10  shows  these  regions  for  £  >  0  and  several  values  of  K 2  >  0.  A  variation  of 
the  wavenumber  from  k  to  k  + 1  results  in  a  very  small  increment  of  K 
so,  in  first  approximation,  we  allow  K2  to  vary  continuously  in  (53)- (54)  and  then 

7f 

the  CW  with  negative  9  {9  e]  —  — ,  0[  in  Fig.  10)  are  all  rendered  unstable.  In  exactly 
the  same  way,  if  £  <  0  then  all  the  CW  with  positive  9  are  destabilized,  see  Fig.  11. 


Figure  10:  Dispersive  instability  conditions  (53)-(54)  for  £  >  0,  a  =  1,  and  different  values  of  K 2 
(arrows  point  towards  increasing  K2  direction  and  shading  indicates  instability). 

Hence,  the  CW  stability  diagrams  in  Figs.  8  and  9  have  to  be  complemented 
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2  4  4  2 

Figure  11:  Dispersive  instability  conditions  (53)-(54)  for  e  <  0,  a  =  5,  and  different  values  of  I\ 2 
(arrows  point  towards  increasing  K2  direction  and  shading  indicates  instability). 

with  the  following  dispersive  instability  criterion: 

-  for  £  >  0  (e  <  0),  all  CW  with  9  <  0  (6  >  0)  are  unstable.  (55) 

In  other  words,  the  dispersive  instability  destabilizes  all  the  CW  along  the  lower 
(upper)  branch  in  Fig.  5  if  £  is  positive  (negative).  And  thus,  for  a  given  power 
inside  the  fiber  p,  there  can  be  bistability  only  for  p  >  pc,  k  <  kc  and  £  <  0  (see 
Fig.  9). 

The  small  dispersion  terms  have  a  dramatic  effect  in  the  stability  of  the  CW  and 
they  should  not  be  ignored.  The  dispersive  instability  is  not  a  higher  order  effect 
that  appears  in  a  much  larger  time  scale;  its  growth  rate  is  of  order  unity  (and 
remains  of  order  unity  as  £  — >  0,  see  eqs.  (51)-(52))  and  it  develops  in  the  transport 
time  scale  ( t  ~  1)  destabilizing  small  scales  with  typical  size  \f\e\  -C  1.  Notice  also 
that,  in  the  limit  £  — >  0,  the  dispersive  instability  condition  (55)  depends  only  on 
the  sign  of  the  dispersion  coefficient  (relative  to  the  nonlinear  coefficients)  and  it 
is  present  for  both  signs.  For  the  CW  with  power  higher  than  p^,  the  beginning 
of  the  destabilization  of  the  intermediate  scales  can  be  seen  from  the  k  — »  oo  limit 
of  the  stability  analysis  without  dispersion  (see  Fig.  7),  but  for  all  other  CW  this 
instability  is  not  detected  if  the  dispersive  terms  are  not  taken  into  account,  and  the 
use  of  the  NLCME  eqs.  (36)-(37)  would  therefore  lead  us  to  wrong  conclusions. 
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5.3  Numerical  simulations  of  the  amplitude  equations 


In  the  following,  we  perform  several  numerical  integrations  of  the  NLCME  equations 
with  small  dispersive  terms  (33)- (35)  in  order  to  check  the  CW  stability  boundaries 
obtained  in  the  previous  sections.  The  details  of  the  numerical  method  can  be  found 
in  Appendix  A. 


Figure  12:  Stability  properties  of  the  CW  (dispersive  instabilities  not  shown)  for  a  =  |  and  k  =  1, 
shading  indicates  instability  and  dots  correspond  to  the  cases  1  and  2  CW. 

We  begin  with  the  case  1  CW  that  corresponds  to  the  parameter  values  o  = 
k  =  1,  p2  =  1  and  9  —  —  |  and,  according  to  eqs.  (41), (42),  has  wavenumber  m  —  0 
and  frequency  a  =  —0.25.  The  reconstructed  original  physical  fields  (see  eq.  (18)) 
produced  by  this  spatially  uniform  CW  take  the  form  of  a  standing  wave  pattern 
along  the  fiber.  As  it  can  be  seen  in  Fig.  12,  the  only  source  of  instability  for  this 
CW  is  the  dispersive  instability:  if  £  <  0  then  the  perturbations  containing  small 
dispersive  scales  are  stable  and  unstable  otherwise.  This  is  corroborated  by  the 
results  of  the  numerical  simulations  of  the  amplitude  equations  (33)- (35)  presented 
in  Figs.  13  and  14,  where  we  plot  the  time  evolution  of  the  norm  of  the  amplitudes, 


for  negative  and  positive  dispersion,  respectively,  starting  from  the  case  1  CW  with 
a  small  perturbation.  For  negative  dispersion  the  case  1  CW  is  stable  and  the  small 
initial  perturbations  do  not  grow  (Fig.  13).  The  perturbations  cannot  decay  to 
zero  because  of  the  absence  of  dissipation  in  the  system,  and  thus  the  value  of  the 
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spatial  derivatives  remains  finite  but  small  (comparable  to  the  size  of  the  initial 
perturbation,  see  Fig.  13).  In  contrast,  if  the  dispersion  is  positive,  the  case  1  CW 
is  unstable  due  to  the  exponential  amplification  of  the  small  dispersive  scales  (see 
Fig.  14).  This  growth  is  not  a  slow  time  effect,  it  takes  place  in  the  time  scale 
t  ~  1  and  remains  constant  as  e  — ►  0,  as  it  can  be  seen  in  the  lower  plot  in  Fig.  14, 
where  we  have  added  for  comparison  the  time  evolution  of  the  spatial  derivatives  for 
two  smaller  dispersion  values.  The  solution  that  develops  consists  of  two  counter- 
propagating  wavetrains  with  dispersive  wavelength  (~  yfe)  and  is  shown  in  Fig.  15; 
notice  how  the  number  of  peaks  is  approximately  doubled  when  the  dispersion  is 
divided  by  4. 
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Figure  13:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (d_)  and  its 
spatial  derivative,  for  a  =  t,  k  =  1  and  dispersion  e  =  — 10-3.  The  initial  condition  is  the  case  1 
CW  with  a  random  perturbation  of  size  ~  10-3. 

In  order  to  check  the  zero  wavenumber  stability  boundary  (47)  we  have  numer¬ 
ically  integrated  the  amplitude  equations  (33)- (35)  starting  from  the  case  2  CW 
(cr  =  |,  k  =  1,  p2  =  4.5,  6  =  —  |  ,  m  =  0  and  a  =  2.375,  see  Fig.  12)  with  a  small 
random  perturbation  and  negative  dispersion,  i.e.,  without  dispersive  instabilities. 
Fig.  16  shows  the  resulting  solution.  The  small  values  of  the  spatial  derivatives  indi¬ 
cate  that,  as  predicted  by  the  linear  stability  theory,  the  destabilizing  perturbation 
mode  is  the  uniform  one  (k  —  0).  Moreover,  the  solution  appears  to  remain  almost 
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Figure  14:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and  its 
spatial  derivative,  for  cr  =  k  =  1  and  dispersion  e  =  10-3  (dashed  lines  correspond  to  smaller 
dispersion  values).  The  initial  condition  is  the  case  1  CW  with  a  small  random  perturbation  of 
size  ~  10-3. 


Figure  15:  Snapshots  of  the  solution  of  the  amplitude  equations  (33)-(35)  at  t  =  60  (|A+|:thick  line, 
|A_|:thin  line),  with  parameters  as  in  Fig.  14  and  dispersion  £  =  1CP3  (above)  and  e  =  10-3/4 
(below). 
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uniform  for  all  times  and  oscillates  periodically  between  two  states:  one  with  both 
amplitudes  nearly  equal  (close  to  a  SW  in  the  original  physical  variables)  and  other 
with  A+  dominating  over  A+  (close  to  a  TW).  Note  that  eqs.  (33)- (35)  preserve  the 
total  energy  of  the  system,  i.e., 

jJ\\A*f  +  \A-\2)dx  =  0, 

and  therefore  if  the  norm  of  one  of  the  amplitudes  increases  the  other  has  to  decrease. 


io~2 


pjll 

10-3 


10->l - 1 - 1 - 1 - 1 - ' - ' - ' - ' - ' - 1 

0  10  20  30  40  50  60  70  80  90  ,  100 


Figure  16:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (T_)  and  its 
spatial  derivative,  for  cr  =  i,  k  =  1  and  dispersion  e  =  — 10-3.  The  initial  condition  is  the  case  2 
CW  with  a  random  perturbation  of  size  ~  10-3. 

In  the  case  3  CW  (a  —  |,  k  —  2,  p2  —  7.2,  6  =  —  f  ,  rri  —  0  and  a  =  3.4),  if  the 
dispersion  coefficient  is  kept  negative,  the  destabilization  is  clue  only  to  condition 
(50),  see  Fig.  17.  The  time  evolution  of  the  solution  of  the  amplitude  eqs.  (33)- 
(35)  starting  from  the  perturbed  case  3  CW  with  £  =  — 10-3  is  shown  in  Fig. 
18.  The  norms  of  the  amplitudes  are  almost  constant  during  the  integration  time, 
but  the  growth  of  the  derivatives  clearly  indicates  the  onset  of  the  instability  that 
takes  the  system  away  from  the  starting  uniform  CW  solution.  The  perturbation 
mode  that  destabilizes  the  CW  is  that  with  k  —  1,  i.e.,  that  exhibiting  one  single 
wavelength  inside  the  domain  (see  Fig.  19),  in  agreement  with  the  linear  stability 
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Figure  17:  Stability  properties  of  the  CW  (dispersive  instabilities  not  shown)  for  a  =  \  and  k  =  2, 
shading  indicates  instability  and  dots  correspond  to  the  cases  3,  4  and  5  CW. 

theory  predictions. 
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Figure  18:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (4“)  and  its 
spatial  derivative,  for  cr  =  i,  k  =  2  and  dispersion  e  =  — 10~3.  The  initial  condition  is  the  case  3 
CW  with  a  random  perturbation  of  size  ~  ICC3. 

For  the  case  4  CW  (a  =  k  =  2,  p2  =  2,  6  =  |  ,  m  =  0  and  a  =  3.5,  see 
Fig.  17)  the  situation  is  the  opposite  of  that  of  the  case  1  CW.  Now  the  CW  is 
stable  for  positive  values  of  the  dispersion,  as  it  can  be  observed  in  the  numerical 
simulations  represented  in  Fig.  20,  and  the  dispersive  instability  develops,  producing 
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Figure  19:  Spatial  profiles  of  the  solution  of  the  amplitude  equations  (33)-(35)  at  times  t  =  0,  5, 6 
and  7,  with  all  parameters  as  in  Fig.  18. 


Figure  20:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and  its 
spatial  derivative,  for  a  =  k  =  2  and  dispersion  £  =  10~3.  The  initial  condition  is  the  case  4 
CW  with  a  random  perturbation  of  size  ~  10~3. 
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Figure  21:  Top:  thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and  its 
spatial  derivative,  for  a  =  -i,  k  =  2  and  dispersion  £  =  — 10~3  (dashed  lines  correspond  to  smaller 
dispersion  values).  The  initial  condition  is  the  case  4  CW  with  a  small  random  perturbation  of 
size  ~  10-3.  Bottom:  space-time  representation  of  the  solution  for  two  time  units  after  t  =  80. 
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Figure  22:  Thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and  its 
spatial  derivative,  for  a  =  k  =  2  and  dispersion  e  =  10~3.  The  initial  condition  is  the  case  5 
CW  with  a  random  perturbation  of  size  ~  10~3. 

small  dispersive  scales  all  over  the  domain,  for  negative  dispersion  (see  Fig.  21).  The 
time  plot  of  ||A±||  in  Fig.  21  shows  the  development  of  the  dispersive  instability 
for  several  dispersion  values;  note  that,  as  predicted  by  the  linear  stability  theory, 
the  dispersive  instability  exponent  tends  to  a  nonzero  constant  as  the  dispersion 
coefficient  is  reduced  £  — ►  0.  A  space-time  representation  of  the  solution,  once 
the  dispersive  scales  are  well  developed,  is  also  included  in  Fig.  21.  Notice  that  the 
small  dispersive  scales  are  largely  advected  by  the  group  velocity  for  short  times,  but 
they  also  evolve  on  a  slower  timescale,  f  ~  1,  to  produce  very  complicated  spatio- 
temporal  dynamics.  And  notice  also  that  the  final  state  in  terms  of  the  original 
physical  variables  is  even  more  complicated:  it  is  the  superposition  of  two  faster 
counter-propagating  wavetrains  with  amplitudes  modulated  by  A+  and  A~  (see  18). 

Finally,  the  case  5  CW  (cr  =  n  =  2,  p2  =  2,  6  =  0.1657  ,  rri  —  1  and 
a  =  7.645)  is  in  the  unstable  side  of  condition  (48),  see  Fig.  17.  The  numerical 
simulation  shown  in  Fig.  22  indicates  that,  again  in  accordance  with  the  linear 
stability  results,  the  CW  is  unstable  and  the  solution  develops  small  dispersive 
scales.  Notice  that  the  growth  rate  of  the  perturbations  is  much  larger  for  this 
case  than  for  the  previous  ones.  This  is  due  to  the  fact  that  the  extension  of  the 
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interval  of  unstable  wavenumbers  is  now  much  longer:  it  starts  in  the  k  rsj  1  (48) 
region  and  extends  up  to  the  dispersive  region  k  i/vFI  (see  Figs.  10  and  11). 
Note  also  that  the  CW  is  inside  the  instability  region  that  corresponds  to  the  +  sign 
in  (48)  and  hence  only  the  amplitude  A+  is  unstable,  as  it  can  be  seen  in  the  initial 
growth  of  the  derivative  in  the  lower  plot  of  Fig.  22. 


37 


6  Gap  Solitons 

Gap  Solitons  (GS)  are  solitary  wave  solutions  of  the  NLCME  (36)-(37)  with  tem¬ 
poral  frequencies  inside  the  frequency  gap  that  opens  up  in  the  linear  propagation 
characteristics  of  the  grating,  where  uniform  infinitesimal  travelling  wave  states  are 
not  possible  (see  Fig.  4).  These  pulse-like  solutions  are  the  result  of  the  combined 
effect  of  nonlinearity  and  grating  reflection  and  have  the  property  that  they  can 
propagate  along  the  grating  with  any  velocity  between  zero  and  the  speed  of  light 
in  the  bare  fibre.  This  is  a  very  promising  feature  that  gives  the  possibility  (at  least 
theoretically)  to  have  localized  packets  of  electromagnetic  energy  at  zero  or  very  low 
velocity  in  the  laboratory  frame. 

This  family  of  solutions,  which  are  not  true  solitons  since  the  NLCME  with  cr  /  0 
are  not  an  integrable  system,  were  introduced  in  [10,  11]  and  can  written  in  the  form 
(see  e.g.  [13,  14,  6]) 

A+s(X,t)  =  -y/KH(£)  e-»/2+vp(0-«7t  (56) 

AqS(x,  t)  =  ^H(i)  ey/2+i^~iKlt  (57) 

where  £  =  k(x  —  ct )  is  a  moving  spatial  coordinate,  c  is  the  velocity  of  the  GS 

c  =  tanhy,  with  —  oo  <  y  <  +oo,  (58) 

^7  is  its  temporal  frequency,  with 

7  =  cos 0/  coshy  =  cos Q\J  1  —  c2,  with  0  <  0  <  tt,  (59) 


and  the  functions  H  and  cp  are  given  by 


m 


2a  sinh(2y)  arc  tan  (tan  (0/2)  tanh[(cosh  y  sin  0)£]} 
1  +  a  cosh  2 y 
sin  0 

a/1  +  a  cosh  2 y  cosh[(cosh  y  sin  0)£  +  i6/ 2] 


+  (sinh  y  cos  0)/ 


(60) 

(61) 


Note  that 

ff(-0  =  H(i),  vH)  =  -¥>(«),  (62) 


and  that  \H\  reaches  its  maximum  at  ^  =  0  and  then  decays  to  zero  monotonically 
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and  exponentially  fast  as  £  — >  ±oo. 

The  GS  are  solutions  of  the  NLCME  defined  for  —  oo  <  x  <  oo  that  correspond 
to  localized  pulses  propagating  with  constant  velocity  c  over  a  zero  background  and 
oscillating  with  frequency  K'y,  like  the  one  represented  in  Fig.  23. 


Figure  23:  Left:  Gap  Soliton  at  t  =  0,  A+{A~ )  thick  (thin)  line,  with  parameters  a  =  k  =  1, 
c  =  i  and  9  =  2.  Right:  domain  of  existence  of  the  GS  in  the  c  —  7  plane  (the  dot  corresponds  to 
the  GS  shown). 

For  fixed  values  of  a  and  k,  this  is  a  two  parameter  (y,  6)  family  of  solutions  with 
velocity  and  frequency  in  the  ranges  —  1  <  c  <  1  and  —y/1  —  c2  <  7  <  \/l  —  c2,  see 
eqs.  (58)  and  (59).  In  the  following  we  will  represent  the  GS  family  in  the  c  —  7 
plane  where  it  fills  an  open  circle  of  radius  one  (Fig.  23).  Notice  also  that  the  GS 
remain  invariant  under  the  changes  c  — >  — c,  x  — >  —x  and  AqS  AqS  and  therefore 
we  can  limit  our  analysis  to  the  GS  with  c  >  0. 
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7  Stability  of  the  Gap  Solitons 

In  this  section  we  compute  the  linear  stability  characteristics  of  the  GS  in  the  con¬ 
text  of  the  NLCME  with  small  dispersion  terms  (33)- (34).  We  will  proceed  as  we 
did  in  Section  5  with  the  CW:  we  will  first  study  the  stability  of  the  GS  against 
perturbations  without  dispersive  scales,  then  the  stability  against  dispersive  pertur¬ 
bations  and  finally  we  will  perform  some  numerical  simulations  to  check  the  stability 
predictions. 

Before  starting  with  the  stability  calculations  it  is  convenient  to  use  the  new  space 
variable  £  =  k{x  —  ct )  to  write  the  amplitude  equations  in  a  frame  moving  with  the 
soliton  and  to  eliminate  the  strength  of  the  grating  k  >  0  and  the  frequency  of  the 
soliton  using  the  new  time  variable  r  =  nt  and  the  new  amplitudes 

71+  =  exp[— iK'yt]. 

With  all  these  changes  the  amplitude  equations  (33)- (34)  take  the  form 

5+  -  (c+  1)5+  =  i£KB+  +i75+  +  iB~  +i(a\B+\2  +  |5~|2)5+,  (63) 

B~  —  (c  —  1)50  =  ienB^  +  i^B~  +  iB+  +  i(<j\B~\2  +  |5+|2)5~,  (64) 

where  |e|  <C  1  and  the  GS  (56), (57)  are  now  steady  states  (up  to  first  order  correc¬ 
tions  in  c) 

Bis  =  -H(0e-y/2+MS\ 

BBs  =  J?(e)e”/2+ivK). 

The  linear  stability  of  the  GS  is  analyzed  as  usual,  taking 

5+  =  5+5  +  b±,  with  |6±|  <  1, 

to  eqs.  (63)-(64)  and  linearizing,  to  obtain  the  system  of  equations 

6+  -  (c  +  1)6+  =  i[e/c6j  +  (7  +  (j)f)b+  +  (1  +  fa)b~  +  0^ b+  +  0 Jr],  (65) 
b~  -  (c  -  1)60  =  i[enb^  +  (7  +  00)6“  +  (1  +  02 )b+  +  03  b~  +  04&+],  (66) 
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with  ^-dependent  coefficients  that  are  given  by 

MO  =  2a|B±s|2  +  IB&I2  =  (e±»  +  2oe™)\H(0\2  (67) 

MO  =  B%SB$S  =  (68) 

<#«)  =  <B*M  =  ae-y*^HM\  (69) 

M(0  =  0Bas?  =  (70) 

MO  =  Bgs'Bqs  =  -e^\H(0\2,  (71) 


and  verify  (see  (62)) 

0f(-£)  =  (0  =  0f  (0,  02(-O  =  02(0, 

03  (“0  =  03(0,  04(-O=04(O- 


7.1  Perturbations  without  dispersive  scales 


In  this  case,  the  perturbations  exhibit  only  scales  that  remain  £  ~  1  as  e  — >  0.  This 
allow  us  to  neglect,  in  first  approximation,  the  dispersion  terms  in  the  system  (65)- 
(66),  namely  we  set  £  =  0  in  eqs.  (65)- (66).  The  resulting  equations  together  with 
their  complex  conjugates  form  a  linear  system  of  four  equations  whose  solutions  can 


be  written  in  the  form 

b+ 1  r  x+(o ' 

5-  =  X-(0 

l+  y+(0 

b-  \  [  Y-( 0 


(72) 


with  uj  and  (X±,Y±)  given  by  the  following  ODE  eigenproblem 


uX+  +  i(c+  1)X+  =  (7  +  0+)X+  +  (l  +  02)X-  +0+r+  +  04y'-,  (73) 

uX~+i{c*  1)X~  =  (7  +  <fq)x~  +  (1  +  02)X+  +  03  Y"  +  04y+,  (74) 

-uY+  -i(c+l)Y+  =  (7  +  0+)F+  +  (1  +  02)Y-  +  0+X+  +  faX~,  (75) 
-uY- -  i(c  -  l)Y~  =  (7  +  07)Y“  +  (1  +  02)F+  +  03  X~  +  04X+.  (76) 

The  spectrum  of  the  problem  above  is  composed  of  all  u6C  with  associated  nonzero 
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eigenfunctions  (X±,Y±)  that  remain  bounded  as  £  — »  ±oo,  either  localized  (point 
spectrum)  or  not  (essential  spectrum). 

A  GS  is  then  unstable  if  there  are  values  of  u  in  its  spectrum  with  negative 
imaginary  part.  Note  that  the  symmetries  of  eqs.  (73)-(76) 

(u,^,X±,Y±)  -  (-u;,-^,^)  (77) 

and 

u  — >  —u),  X±  Y±  (78) 

allow  us  to  restrict  the  analysis  to 

>  0  and  Qw  >  0. 

The  first  symmetry  above  follows  from  the  spatial  reflection  (23)  and  time  reversal 
(30)  invariance  of  the  system  and  the  second  is  due  to  the  particular  form  of  the 
perturbation  (72). 

In  the  limit  |£|  — >  oo  the  system  (73)-(76)  becomes  a  constant  coefficient  linear 
system  and  its  solutions  can  be  written  in  the  form 

(X±,Y±)  =  (Vjf.iyje**  (79) 

where  the  complex  constants  X^t  Y(^  and  a  satisfy 

[y  —  o’  +  (c  +  l)a]A0h  +  Ag  =  0,  (80) 

[7  —  u  +  (c  —  1)ck]  Ag  +  Xq  =  0,  (81) 

and 

[7  +  u  +  (c  +  l)a]Y0+  +  Y0  =  0,  (82) 

[7  +  cu  +  (c  —  l)a]Vg  +  Tg+  =  0.  (83) 

Imposing  that  these  systems  exhibit  nontrivial  solutions,  we  obtain 

(7  —  u  +  ca )2  —  a2  —  1  =  0  (84) 
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for  (80)-(81)  and 

(7  +  ijj  +  ca )2  —  a2  —  1  =  0  (85) 

for  (82)-(83),  and  then  for  every  value  of  c 0  we  have,  in  general,  four  complex  values 
of  a  that  give  the  behavior  of  the  eigenfunctions  as  £  — >  ±00. 

Real  values  of  a  correspond  to  eigenfunctions  with  oscillatory  tails  that  remain 
bounded  and  do  not  decay  to  zero  as  £  — >  ±00,  see  eq.  (79).  From  eqs.  (84)  and  (85), 
real  values  of  a  are  seen  to  occur  only  for  real  values  of  to  along  the  four  branches 

u  =  ca  ±  7  ±  Va2  +  1,  afl. 

This  defines  the  essential  spectrum  of  (73)-(76),  which  is  confined  to  the  real  line  and 
therefore  it  does  not  give  rise  to  any  instabilities.  The  essential  spectrum  is  composed 
of  four  semi-infinite  segments  of  the  real  line  bounded  by  the  four  extrema  of  the 
expressions  above:  ±uo  and  (see  Fig.  24)  ,  where  uq  >  0  and  oj\  >  0  are  given 
by 

u0  =  —7  +  Vl  —  c2  and  uq  =  7  +  a/1  —  c2. 


1 0 

/-\l\  +1/1  -c2 

/ 

\  |7|+/T^2 

Figure  24:  Essential  spectrum  in  the  complex  to  plane  (thick  line)  with  dots  marking  its  edges. 


Away  from  the  essential  spectrum,  the  roots  of  eqs.  (84)  and  (85), 

~c{ u  ~  7)  ±  a/(cu  -  y)2  -  (1  -  c2) 


and 


ax  — 


c^y  — 


1  —  c2 

— c( U)  +  7)  ±  y/(u  +  7)2  -  (1  -  C2) 

1  —  C2 


(86) 


(87) 


respectively,  have  nonzero  imaginary  parts  with  opposite  signs  (the  +  (— )  super 
index  corresponds  to  the  root  with  positive  (negative)  imaginary  part).  This  is 
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readily  seen  after  taking  into  account  that  the  imaginary  parts  of  the  roots  a \  can 
be  regarded  as  two  real  functions  that  depend  continuously  on  ou,  do  not  vanish  in 
the  connected  set  obtained  after  removing  the  essential  spectrum  from  the  complex 
plane  (recall  that  the  essential  spectrum  is  precisely  the  set  where  Tf a  =  0)  and  have 
different  signs  at  u  —  0;  and  the  same  holds  true  for  a?.  In  order  to  ensure  that 
both  |A"±|  and  F±|  remain  bounded,  we  must  remove  the  exponentially  growing 
eigenfunctions,  and  the  only  allowed  asymptotic  behaviors  as  £  — >  ±oo  are 

(X+,  X-)  =  (-1, 7  -  u;  +  (c  +  1  )a±)eia*?  (88) 

(F+,  Y~)  =  (-1, 7  +  u  -  (c  +  1  )a$)eia^.  (89) 

Summarizing,  outside  the  essential  spectrum,  the  eigenfunctions  are  given  by 
(73)-(76)  together  with  the  boundary  conditions  (88)  and  (89). 

The  point  spectrum  is  now  calculated  by  means  of  the  so-called  Evans  function, 
which  is  defined  as  follows  (see  e.g.  [32,  33]  and  [34]  for  a  recent  review  on  this 
topic).  For  each  w  6  C  outside  the  essential  spectrum  we  integrate  (73)-(76)  both 
from  £  — >  — cx)  and  from  £  — >  +oo  and  starting  from  a  linear  combination  of  (88) 
and  (89),  to  obtain 

(Xt  y-)  =  MX?,  Y±)  +  k2(x?  ,  y±) 

and 

(At,  Ft)  =  A'3(  At,  It)  +  At  At,  Y?), 

respectively.  Note  that  once  the  appropriate  asymptotic  behavior  is  imposed,  each 
orbit  depends  linearly  on  two  (complex)  constants.  These  two  orbits  must  have  a 
non  void  intersection  at  £  =  0,  namely  they  must  satisfy  (X?,Yj?)  =  (X?,Y+) 
at  £  =  0.  This  is  a  homogeneous,  fourth  order  linear  system  for  K\, ,  K4  and 
we  define  the  Evans  function,  E(u),  as  the  value  of  the  determinant  its  coefficient 
matrix.  Therefore,  it  has  nontrivial  solutions  if  and  only  if 

E(u)  =  0. 

In  other  words,  the  Evans  function  detects  the  intersections  of  the  stable  and  unsta- 
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ble  manifold  of  the  fixed  point  (X^,  Y^  =  (0,0)  of  (73)-(76)  and  it  is  defined  only 
when  (0,  0)  is  hyperbolic,  that  is,  for  u  values  outside  the  essential  spectrum.  It  is 
an  analytic  function,  as  it  follows  from  the  analyticity  of  the  constituent  solutions 
(X^,  Y^), . . .  (Xf,  Y^),  with  the  property  that  its  zeros  coincide  with  the  eigenval¬ 
ues.  Note  that  because  of  the  symmetries  (77)  and  (78)  the  Evans  function  must 
also  be  such  that 

E{—u)  =  E(u )  and  E{u)  =  E{u).  (90) 

We  can  now  analyze  the  stability  of  the  GS  just  by  applying  the  principle  of 
the  argument  to  the  Evans  function  along  the  contour  depicted  in  Fig.  25,  which 
encircles  the  half  complex  plane  that  corresponds  to  c o  with  positive  imaginary  part, 
to  obtain  the  number  of  unstable  eigenvalues. 


Figure  25:  Contour  in  the  complex  to  plane  for  the  application  of  the  principle  of  the  argument  to 
E(u).  The  thick  line  correspond  to  the  essential  spectrum  and  the  black  dots  mark  its  edges. 

In  order  to  do  this  we  have  to  take  into  account  the  following  considerations: 

•  Using  the  analytic  extension  of  the  square  roots  in  (86)  and  (87),  we  can  extend 
analytically  the  Evans  function  from  the  half  complex  plane  Qw  >  0  to  the 
essential  spectrum,  except  at  its  edges,  ±uq  and  ±tui,  where  the  extended 
Evans  function  presents  singularities.  Note  that  the  Evans  function  is  not 
continuous  at  the  essential  spectrum  and  thus  its  extensions  from  above  and 
below  (namely,  from  ‘Acu  >  0  and  from  Qtn  <  0)  do  not  coincide. 

•  The  number  of  unstable  eigenvalues  is  given  by 

n  =  ^-A  fa,  (91) 
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where  AiJ>c  is  the  increment  of  the  argument  of  E  along  the  contour  in  Fig.  25 
(see  e.g.  [35]  for  a  general  description  of  the  principle  of  the  argument),  which, 
due  to  the  first  symmetry  in  (90),  can  be  written  as  twice  the  increment  of  the 
argument  along  the  right  half  of  the  contour. 

•  u  =  0  is  an  eigenvalue  with  geometric  multiplicity  2  and  associated  eigenfunc¬ 
tions 


(V±,r±)  =  (B*sf,B±S{), 

(v±,y±)  =  (iB±s,-!B±s), 

which  result  from  the  two  continuous  symmetries  of  the  system  (63)- (64) 

+  and  BA  — >  eiC2B±.  (92) 

But  its  algebraic  multiplicity  is  four  because  we  have  also  two  generalized  eigen¬ 
functions 


(^±,  Y±)  =  {-idB±s/dc,idB±s/dc), 

(X±,Y±)  =  (-idB^/d^idB^/dj), 

which  come  from  the  fact  that  b±  =  d B^g/dc  and  =  dB^g/d 7  satisfy 

-Best  ~  (c  +  l)6f  =  iii'y  +  (pt)b+  +  (1  +  </>2)&_  +  <j>tb+  +  <f>Jr], 
-BGsa  -  (c-  1)6^  =  *[(7  +  4>i )b~  +  (1  +  02)&+  +  03^  +  <M+], 

and 


-iBcs  ~  (c  +  1  )6f  =  *[(7  +  (j)f)b+  +  (1  +  ( fo)b  +  ^ b+  +  fab  ], 

-iBGS  -  (c-  1  )bf  =  i[( 7  +  4>i)b~  +  (1  +  (f)2)b+  +  4>^b~  +  04^+], 

respectively,  as  it  can  be  seen  upon  differentiation  with  respect  to  the  soliton 
parameters  c  and  7  in  the  GS  equation,  i.e.,  in  the  steady  state  version  of  (63)- 
(64)  with  e  =  0.  In  other  words,  the  algebraic  multiplicity  is  equal  to  four  due  to 
the  four-dimensionality  of  the  family  of  gap  solitons  (that  can  be  parametrized 
in  terms  of  c,  7  and  the  two  parameters  in  the  continuous  symmetries  (92)), 
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but  c  and  7  do  not  produce  any  increment  of  the  geometric  multiplicity  because 
of  the  time-dependent  character  of  change  of  variables  performed  in  (63)- (64). 
The  contribution  to  (91)  of  the  fourth  order  root  u  =  0  is  thus  equal  to  —2,  see 
Fig.  25. 

•  The  evaluation  of  the  contribution  of  the  outer  arc  in  Fig.  25  requires  to  know 
the  asymptotic  behavior  of  E  as  |cu|  — >  00.  In  this  limit,  eqs.  (73)-(76)  simplify 
drastically  and  the  allowed  solutions  according  to  (88)- (89)  are  given,  at  first 
order,  by 


(_e^/(i+c),  o,o,0),  (0,0,— eiaJ«/(1+c),0)  for  £>0  and 

(0,  e/{1"c),0,0),  (0,0,0,  ^e-^/M)  for  £  <  0, 

recall  that  Thu  >  0.  Therefore,  E(u)  ~  cu2  as  |u;|  — >  00,  and  the  arc  gives  a 
contribution  of  amount  1  to  (91). 

•  The  second  symmetry  in  (90)  ensures  that  the  Evans  function  is  real  in  the 
segment  of  the  real  line  that  is  not  part  of  the  essential  spectrum  (labeled  0  in 
Fig.  25)  and  its  contribution  to  (91)  is  equal  to  —n0/ 2,  where  n0  is  the  number 
of  roots  along  the  segment. 

•  The  extended  Evans  function  is  not  analytic  at  the  edges  of  the  essential  spec¬ 
trum  but  typically  they  do  not  contribute  to  (91)  because  E  ^  0  at  u>  —  ±cuo,i 
in  a  generic  situation.  The  same  argument  can  be  used  to  conclude  that,  gener- 
ically,  there  are  no  zeros  in  the  rest  of  the  essential  spectrum  (segments  1  and 
2  in  Fig.  25). 

Collecting  all  the  above  mentioned  contributions,  the  number  of  unstable  eigenvalues 
can  be  finally  written  as 

n  =  —n0  -  1  +  -(A 1/7  +  A-02),  (93) 

7 r 

where  A^1j2  represent  the  increment  of  the  argument  along  the  segments  1  and  2  in 
Fig.  25.  Note  that  the  formula  above  requires  the  evaluation  of  the  Evans  function 
only  along  the  positive  real  axis. 

The  stability  of  the  GS  has  been  numerically  analyzed  using  (93)  for  different 
values  of  the  parameters  c  and  7.  The  evaluations  of  the  Evans  function  have  been 
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performed  with  an  adaptive  Runge-Kutta  method  of  4th  order  and  the  integrations 
have  been  started  at  £  =  ±£00,  where  £oo  has  been  taken  large  enough  to  ensure  that 
the  GS  amplitude  is  negligible  at  this  point. 

The  results  we  have  found  are  represented  in  Fig.  26  and  can  be  summarized 
as  follows:  (i)  the  GS  are  unstable  for  values  of  7  smaller  than  a  certain  7C  that 
is  negative  very  small  (see  Fig.  26),  and  (ii)  the  onset  of  the  instability  always 
takes  place  through  the  collision  of  two  real  eigenvalues  in  the  0  segment  in  Fig.  25 
that  become  complex  with  a  nonzero  value  of  (temporal  frequency)  that  is  also 
presented  in  Fig.  26.  For  7  <  yc  more  eigenvalues  become  unstable,  both  with  zero 
and  nonzero  temporal  frequencies,  but  the  stability  is  never  regained. 


Figure  26:  Left:  GS  stability  diagram  (shading  indicates  instability  and  dots  correspond  to  the 
numerical  simulations  in  Section  7.3).  Right:  frequency  of  the  neutral  mode  at  instability  onset. 

These  stability  results  agree  with  those  presented  by  Barashenkov  and  collabo¬ 
rators  in  [13,  14],  which  were  calculated  using  Newton’s  method  combined  with  a 
finite  difference  approximation  on  a  finite  interval  of  eqs.  (73)- (76). 

7.2  Perturbations  with  dispersive  scales 

We  now  complete  the  stability  analysis  of  the  GS  with  the  study  of  the  effect  of 
perturbations  that  contain  small  dispersive  scales,  that  is,  scales  with  a  typical  size 
that  goes  as  y/jej  when  e  — >  0. 

In  this  case,  the  dispersion  terms  cannot  be  neglected  in  eqs.  (65)- (66)  and  the 
eigenvalue  problem  that  is  obtained,  after  assuming  exponential  time  dependence 
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(see  eq.  (72)),  can  be  written  as 

wl+  +  i(c  +  1)A+  =  £KX+  +  (7  +  <f>t)X++  (1  +  <f>2)X-+  0+A++  0 4y-,  (94) 
uX~  +  i(c  -  1)A-  =  £kX~  +  (7  +  07)X-+  (1  +  02)A++  faY~+  04A+,  (95) 
-cuF+-  i(c  +  1)A?+  =  £/7T+  +  (7  +  <t>t)Y++  (1  +  02 )Y~+  $tx++  04*“,  (96) 
-c oY~  -  i{c  -  1  )Yf  =  ekY#  +  (7  +  07 )Y~+  (1  +  02)y++  0JA-+  04A+,  (97) 

where  the  eigenfunctions  [X±,Y±)  must  decay  to  zero  as  £  — >  ±00. 

If  we  now  introduce  the  fast  (dispersive)  spatial  variable  rj  =  7/  and  expand 
the  solution  of  the  system  above  as 

(x±,y±)  =  (At,  roPo.O  +  TIdtv?,  57)0,  ?)  +  •••, 

ui  =  oj0/ vTi  +  wi  + . . . , 

with  e|  <C  1,  we  obtain,  at  leading  order,  the  following  uncoupled  homogeneous 
system 


<^0*0"  +  i(c  +  1)A^ 

=  0, 

(98) 

uoxo  +  i(c  -  1)A0“ 

=  0, 

(99) 

-u0Y0+  -  i(c  +  1)Y+ 

=  0 

(100) 

-u0Y0~  -  i(c  -  1  )Y0~ 

=  0, 

(101) 

whose  solutions  are  of  the  form 

xo  =  ^0  (6e^/(c±1), 

=  5±(£)eiaw/(c±1), 

where  has  to  be  purely  real  for  the  eigenfunctions  to  be  bounded  in  the  short  scale 
77.  These  expressions  represent  travelling  waves  that  propagate  in  7  to  the  left  (+) 
and  right  (“),  with  velocities  1  +  c  and  1  —  c,  respectively,  and  exhibit  modulations 
in  the  longer  scale  £  rv  2. . 

At  next  order,  a  system  is  obtained  for  (A^,  Y^)  that  is  identical  to  (98)-(101) 
but  with  a  nonzero  right  hand  side.  The  equations  for  (Aq  (£),  I^(£))  result  from 
applying  solvability  conditions  to  this  system,  i.e. ,  from  forcing  (A±,  Y^)  to  remain 
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bounded  in  the  fast  scale,  and  are  given  by 


+  i(c  +  1)Aq£ 

w\Aq  +^(c  —  1)^40^ 
-urf+^iic+l  )B+ 
-cni^o  -  i(c  -  1)B~£ 


(7  +  0+ 

(7  +  0r 
(7  +  0+ 
(7  +  0r 


- — ;K- 

\e\  1 

[c+  l)2 

£ 

— tK- 

\e\  1 

(c-l)2 

£ 

- — -K- 

£  1 

:c+i)2 

£ 

—  K- 

e\  (c  —  l)2 


Mo  +  0+5o+ , 

)A)  +  03  -B0  , 

W  +  03  ■ ^0  , 
)Bq  +  03  A0  , 


(102) 

(103) 

(104) 

(105) 


together  with  the  conditions  that  Aq  — >  0  and  — >  0  as  £  — >  ±cx). 

Notice  that,  as  it  happened  for  the  CW  in  Section  5.2,  dispersive  perturbations 
propagating  to  the  left  (+)  and  right  (“)  are  uncoupled  and  the  system  above  is 
composed  of  two  second  order  independent  subsystems:  (102), (104)  and  (103), (105). 

In  the  limit  £  — >  ±00,  we  have  0+  — >  0  and  £>3  — >  0  (see  (67)  and  (69))  and 
the  first  subsystem  becomes  a  constant  coefficient,  diagonal  system  with  associated 
spatial  exponents  A, 

Oy,Bo+)  =  (a,/3)e’A£, 


of  the  form 

a  =  (Wl  ±  (7  -  + 

|£|  (c+lffi 

Both  spatial  exponents  have  the  same  imaginary  part  and  for  Qwi  7^  0  the  two 
corresponding  asymptotic  behaviors  grow  simultaneously  unbounded  as  either  £  — ► 
+00  or  £  — >  — cx)  (depending  on  the  sign  of  Qwi),  and  therefore  there  is  no  other 
possible  bounded  solution  apart  from  the  zero  solution.  The  same  holds  true  for  the 
second  subsystem  and  allows  us  to  conclude  that  there  are  no  instabilities  associated 
with  the  perturbations  with  small  dispersive  scales. 

In  contrast  to  what  was  obtained  for  the  CW  in  Section  5.2,  the  GS  do  not  exhibit 
dispersive  instabilities.  The  required  spatial  decay  of  the  dispersive  perturbations  as 
£  — >  ±00  can  only  be  achieved  through  the  coupling  that  comes  from  the  nonlinear 
interaction  terms  in  the  NLCME  (see  eqs.  (102)-(105)),  which  remains  nonzero  for 
the  CW  but  vanishes  in  the  GS  case  as  £  — »  ±00. 
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7.3  Numerical  simulations  of  the  amplitude  equations 


The  GS  stability  predictions  are  cross-checked  here  against  the  results  of  the  numer¬ 
ical  simulation  of  the  NLCME  with  small  dispersion  terms  (33)- (35). 

The  code  developed  for  the  numerical  integration  of  equations  (33)- (35)  assumes 
periodic  boundary  conditions  (see  Appendix  A).  In  order  to  simulate  the  NLCME  in 
an  infinite  domain  we  take  a  spatial  period  large  as  compared  with  the  characteristic 
length  of  the  GS  (Loo)  and  we  then  rescale  space,  time  and  the  size  of  amplitudes 
to  make  the  spatial  period  one  and  the  group  velocity  and  the  first  nonlinear  coef¬ 
ficient  also  equal  to  one  (see  Fig,  27  for  the  explicit  expressions  of  the  rescalings). 
Note  that  after  these  changes  of  variables  the  resulting  strength  of  the  grating,  k, 
is  increased  and  the  dispersion,  e,  is  reduced.  Note  also  that  the  validity  of  the  nu¬ 
merical  simulations  is  limited  in  time  because  the  finite  length  effects  will  eventually 
contaminate  the  solution  on  the  whole  domain. 


K-oo  b  £oc;  l^ool 


Xqq  0  1 


Figure  27:  Rescalings  performed  to  transform  the  large  finite  domain  of  size  L0 0  into  [0, 1]. 

We  present  the  results  of  four  simulations  with  parameters  a  =  k  —  40  and  £  = 
0  (dispersion  is  not  considered  in  this  first  set  of  simulations  because  we  are  interested 
only  in  instabilities  due  to  perturbations  without  dispersive  scales)  starting  from  the 
following  GS  with  a  small  random  perturbation  of  size  ~  10“3  superimposed: 

•  GS1:  c  =  0  and  7  =  cos  1  =  .5403  . . .. 

•  GS2:  c  =  0  and  7  =  cos  2  =  —.4165  . . . 

•  GS3:  c  =  0.5  and  7  =  ^  cos  1  =  .4679  . . . 
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.  GS4:  c  =  0.5  and  7  =  ^  cos  2  =  -.3604  . . . 

According  to  the  stability  analysis  of  the  previous  section,  the  GS1  is  a  stable 
soliton  (7  >  0  in  Fig.  26).  The  simulation  results  are  consistent  with  this:  the 
norm  of  the  solution  and  its  derivatives  is  almost  constant  for  50  time  units  with 
only  a  small  variation  due  to  the  initial  perturbation  (see  top  plot  in  Fig.  28),  and 
the  x  —  t  diagram  in  Fig.  28  indicates  that  the  spatial  profiles  of  the  modulus  of 
the  amplitudes  remain  practically  unaltered,  as  it  should  be  for  a  GS  with  zero 
propagation  velocity. 

On  the  other  hand,  the  GS2  is  unstable  (see  Fig.  26)  and  the  rapid  development 
of  the  oscillatory  instability  can  be  clearly  appreciated  from  the  simulation  results 
plotted  in  Fig.  29.  From  the  evolution  of  the  norm  of  the  solution  in  the  top  plot 
of  Fig.  29  the  frequency  of  the  most  unstable  mode  can  be  estimated  ^estimated  — 
1.34  ...  ,  and  this  number  is  in  very  good  agreement  with  the  real  part  of  the  unstable 
eigenvalue,  oj  ~  1.33  . . .  +  i0.0338  . . .,  computed  applying  Newton’s  method  to  find 
the  zero  of  the  Evans  function.  This  is  actually  the  only  unstable  eigenvalue  as  it 
can  be  seen  in  the  map  of  the  zeros  of  the  Evans  function  that  is  also  included  in 
the  middle  plot  of  Fig.  29. 

The  GS3  is  a  stable  soliton  (see  Fig.  26)  that  now  propagates  with  velocity 
0.5,  and  the  results  of  the  simulation  presented  in  Fig.  30  corroborate  again  the 
conclusions  of  the  stability  analysis:  the  norms  remain  almost  constant  and  the 
x  —  t  diagram  shows  the  soliton  propagating  to  the  right  without  any  perceptible 
distortion  after  50  units  of  time  (note  that  now  the  amplitudes  of  the  GS  are  not 
equal  because  this  is  a  travelling  pulse,  y  ^  0  in  (56)-(58)). 

And  for  the  last  case,  the  GS4,  the  stability  results  indicate  that  it  corresponds 
to  an  unstable  propagative  soliton  (see  Fig.  26)  and  this  is  in  good  agreement  with 
the  simulation  results  presented  in  Fig.  31.  The  instability  that  develops  is  again 
oscillatory  and  its  frequency  estimated  from  the  norm  evolution  of  the  solution, 
^estimated  —  1.17. . .,  is  again  very  close  to  the  unstable  eigenvalue  obtained  solving 
the  equation  E(u>)  =  0,  u>  ~  1.166  . . .  +  i0.0218  . . .  (middle  plot  of  Fig.  31).  Also, 
the  growth  rate  computed  from  E(u)  =  0  is  now  smaller  than  in  the  GS2  case  and 
produces  a  slower  instability  development,  as  it  can  be  appreciated  from  the  time 
plots  of  ||A±||  in  Figs.  29  and  31. 

Finally,  in  order  to  complete  the  numerical  checking  of  the  stability  results,  we 
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have  performed  several  simulations  of  the  stable  solitons  GS1  and  GS3  with  k  =  40 
and  different  values  of  the  dispersion  coefficient,  namely,  £  =  ±4  10”4,  ±2  10-4  and  ± 
10~4.  The  numerics  corroborate  the  theoretical  results  obtained  in  Section  7.2,  that 
is,  no  sign  of  any  dispersive  instabilities  was  detected.  The  numerical  solutions  found 
were  indistinguishable  from  those  represented  in  Figs.  28  and  30  and  therefore  their 
plots  are  omitted. 
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Figure  28:  Top:  thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and 
its  spatial  derivative,  for  a  =  k  =  40  and  £  =  0.  The  initial  condition  is  the  GS1  with  a  random 
perturbation  of  size  ~  10-3.  Bottom:  space-time  representation  of  the  solution  for  four  time  units 
right  after  t  =  50  in  the  top  plot. 
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Figure  29:  Top:  thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and 
its  spatial  derivative,  for  er  =  i,  k  =  40  and  e  =  0.  The  initial  condition  is  the  GS2  with  random 
perturbation  of  size  ~  10-3.  Middle:  $tE(uj)  =  0  (SsE(lj)  =  0)  contours  in  the  complex  oj  plane  in 
solid  (dashed)  line.  Bottom:  space-time  representation  of  the  solution  for  the  last  four  time  units 
in  the  top  plot. 
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Figure  30:  Top:  thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and 
its  spatial  derivative,  for  a  =  k  =  40  and  e  =  0.  The  initial  condition  is  the  GS3  with  random 
perturbation  of  size  ~  10-3.  Bottom:  space-time  representation  of  the  solution  for  four  time  units 
right  after  t=50  in  the  top  plot. 
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Figure  31:  Top:  thick  (thin)  lines  indicate  the  time  evolution  of  the  spatial  norm  of  A+  (A-)  and 
its  spatial  derivative,  for  cr  =  1,  n  =  40  and  e  =  0.  The  initial  condition  is  a  GS4  with  random 
perturbation  of  size  ~  10-3.  Middle:  $tE(uj)  =  0  (SsE(uj)  =  0)  contours  in  the  complex  oj  plane  in 
solid  (dashed)  line.  Bottom:  space-time  representation  of  the  solution  for  the  last  four  time  units 
in  the  top  plot. 
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8  Numerical  simulations  of  the  MLE 


In  order  to  confirm  the  validity  of  the  amplitude  equations  with  small  dispersion 
terms  and  the  theoretical  stability  results  obtained  in  the  previous  sections,  we  have 
tested  them  against  numerical  simulations  of  the  complete  MLE  (16)-(17). 

The  MLE  simulations  correspond  to  a  ring  shaped  fiber  grating  geometry  with  a 
large  number  of  refractive  index  oscillations,  that  is,  to  periodic  boundary  conditions, 


E{x  +  L,t)  =  E(x,  t ),  P(x  +  L,t)  =  P(x,  t ), 


with  I  >  1,  see  Appendix  B  for  the  details  of  the  numerical  method.  The  initial 
conditions  used  in  all  the  simulations  are  slightly  perturbed  CW  and  GS  configura¬ 
tions. 

The  relation  between  the  solutions  of  the  MLE  (16)-(17)  and  the  amplitude  equa¬ 
tions  (33)-(35),  after  taking  into  account  the  rescalings  in  (32),  can  be  written  as 


E 

P 


1  —  cu2 


,  l-^-(A+eix+iut  +  A~e~ix+iut)  +  c.c.  +  . . . .  (106) 

V  L\U2\ 


The  grating  strength  and  the  dispersion  are  related  to  their  scaled  counterparts  by 

An  =  T\9  ,  k  and  £  =  —  (107) 

L\w\  Lvg 

where  L  1  and  the  parameters  u,  vg ,  d,  w  and  «2  are  given  by  (15),  (24),  (25), 
(26)  and  (27)  as  functions  of  cu2  and  only. 

The  structure  of  the  actual  physical  patterns  that  result  from  (106)  becomes  more 
clear  if  we  rewrite  the  electric  held  as 


E  ~  ( A+eiu}t  +  A~e~iu>t)eix  +  c.c.  +  . . . , 

where  the  amplitudes  A±  depend  slowly  on  space  and  time  and  thus  remain  approx¬ 
imately  constant  in  any  region  of  extent  Ax  ~  1  and  for  any  time  interval  At  ~  1. 
The  complex  number  inside  parentheses  in  the  expression  above, 

A+eiu}t  +  A~e~iut  =  Mei<p, 
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can  be  regarded  as  the  sum  of  two  constant  complex  numbers  that  rotate  with 
opposite  angular  velocity  u>  (Fig  32a)  and  its  modulus,  M,  gives  the  local  amplitude 
of  the  envelope  of  the  pattern, 

E  ~  Meiip+ix  +  c.c.  + - M  cos  (a;  +  <p)  +  ..., 

see  Fig  32b,  that  oscillates  in  time  between  the  extreme  values 

Mmax  =  \\A+\  +  \A~\\  and  Mmin  =  ||A+|  -  \A~\\. 

Therefore,  if  the  amplitudes  are  nearly  equal  in  a  certain  region,  jA+|  ~  |A“|,  then 
M  will  oscillate  between  0  and  Mmax  and  ip  will  alternate  between  ±(cpA+  +  <Pa-)/ 2 
(see  Fig32a),  and  the  pattern  will  locally  look  like  a  SW.  On  the  other  hand,  if  one  of 
the  amplitudes  dominates,  |A+|  3>  |v4—  | ,  then  M  is  almost  constant  and  ip  behaves 
approximately  like  cut,  and  the  pattern  resembles  a  TW.  For  any  other  configuration 
with  arbitrary  A+  and  A~  the  resulting  pattern  will  look  more  complicated;  it  will 
be  a  superposition  of  SW  and  TW. 

a)  b) 

E  M  cos(x  +  p) 

A+ 

0 

Figure  32:  Construction  of  the  solution  of  the  MLE  from  the  amplitudes  A±. 

All  CW  simulations  (shown  in  Figs.  33  to  37)  have  been  carried  out  using  the 
following  parameters:  =  1,  Uq  =  2  and  L  =  64(2n)  (i.e.,  there  are  128  oscillations 

of  the  grating  inside  the  fiber  ring).  For  every  simulation,  the  corresponding  value 
of  An  is  calculated  from  k  using  eq.  (107)  and  the  sign  of  the  dispersion  is  selected 
by  either  choosing  the  +  or  the  —  sign  in  the  equation  for  cu  (15),  see  Fig.  3. 
The  number  of  Fourier  modes  and  the  time  step  used  in  the  numerical  scheme  are 
Aif Fourier  =  512  and  At  =  .01  (see  Appendix  B). 
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In  the  results  presented  in  Fig.  33  the  initial  condition  is  given  by  the  CW1  in 
Fig.  12  with  a  small  perturbation  superimposed.  The  first  and  second  plot  show 
the  time  evolution  of  the  norm  of  the  electric  held, 


and  correspond,  respectively,  to  u>  =  u>~  and  u  =  cu+,  that  is,  to  negative  and  positive 
dispersion  (see  eq  (107)  and  Fig.  3).  The  electric  held  pattern  for  the  CW1  is  a 
uniform  SW.  The  associated  value  of  ||i7||  oscillates  very  fast  in  time  between  0  and 
its  maximum  value,  but  this  oscillation  is  too  fast  to  be  appreciated  in  the  hrst  two 
plots  of  Fig.  33  and  they  just  look  like  black  patches.  The  MLE  simulations  confirm 
the  theoretical  predictions:  the  CW1  is  stable  for  negative  dispersion  (hrst  plot  in 
Fig.  33)  and  unstable  for  positive  dispersion  (something  starts  to  develop  around 
t  =  60000  in  the  second  plot).  The  corresponding  spatial  prohles  of  E  at  t  =  75000 
are  given  in  the  third  and  fourth  plots  of  Fig  33:  for  negative  dispersion  (third 
plot)  a  perfectly  uniform  oscillatory  pattern  is  obtained  but,  for  positive  dispersion, 
a  modulation  is  clearly  present  (fourth  plot).  In  order  to  be  sure  that  this  is  a 
dispersive  instability  we  have  repeated  the  unstable  (ui  =  lv+)  MLE  simulation 
in  a  four  times  longer  domain  (L  =  256  (27r)).  The  resulting  spatial  profile  of  E 
at  t  =  160000  is  shown  in  the  last  plot  of  Fig.  33.  Notice  how  the  number  of 
basic  wavelengths  is  now  four  times  higher  but  the  number  of  wavelengths  of  the 
modulation  only  approximately  doubles  (increases  from  5  to  9),  indicating  the  onset 
of  dispersive  scales  whose  size  goes  as  a J~L  ~  4j. 

Dispersion  is  negative  in  the  simulation  shown  in  Fig.  34,  which  starts  from  a 
small  perturbation  of  the  CW2  in  Fig.  12.  According  to  the  stability  results  from 
Section  5,  this  CW  is  unstable  and  evolves  to  a  spatially  uniform,  time  periodic 
solution  (see  Fig.  16).  The  solution  of  the  MLE  reproduces  this  behavior:  the 
fast  oscillations  of  ||Ej|  in  the  top  plot  slowly  alternate  between  going  down  to  zero 
(nearly  equal  amplitudes)  and  remaining  away  from  zero  (different  amplitudes),  and 
E  does  not  develop  any  slow  spatial  modulation  (see  bottom  plot  of  Fig.  34). 

Fig.  35  corresponds  to  a  MLE  simulation  starting  from  the  CW3  in  Fig.  17  and 
with  negative  dispersion.  The  linear  theory  now  predicts  a  destabilization  due  to  a 
mode  with  one  single  wavelength  inside  the  ring.  This  is  again  corroborated  by  the 
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MLE  simulation  results,  as  it  can  be  appreciated  from  the  bottom  plot  of  Fig.  35 
that  shows  the  spatial  profile  of  E  at  the  beginning  of  the  instability  development. 
For  longer  times,  the  ||£'||  plot  in  Fig.  35  indicates  that  the  solution  becomes  more 
complicated. 

The  initial  condition  for  the  two  MLE  simulations  in  Fig.  36  is  the  CW  defined 
by  the  parameters  k  —  1,  9  =  j  and  p2  =  1  (see  Fig.  12)  with  a  small  perturbation. 
The  first  and  third  plot  correspond  to  positive  dispersion  and  indicate  that  the 
CW  is  stable.  The  dispersion  is  negative  in  the  second  and  fourth  plot  where  the 
destabilization  of  the  CW  can  be  clearly  seen  in  the  time  evolution  of  ||F||  and  the 
spatial  profile  of  E  shows  dispersive  modulations.  This  is  in  perfect  agreement  again 
with  the  linear  stability  results  in  Section  5. 

Fig.  37  shows  the  last  of  the  CW  related  MLE  simulations.  Dispersion  is  taken 
positive  in  this  case  and  the  initial  condition  is  a  perturbed  CW  with  parameters 
k  —  2,  6  —  j  and  p 2  =  3,  which  corresponds  to  a  point  in  Fig.  17  above  the 
point  labeled  4  and  inside  the  shaded  unstable  region.  As  the  stability  calculations 
predicted,  the  solution  develops  dispersive  scales  (see  lower  plot  in  Fig.  37)  and  the 
development  of  the  instability  is  now  faster  than  in  the  previous  cases  (top  plot  in 
Fig.  37)  because  now,  as  we  saw  in  Section  5,  the  range  of  unstable  wavenumbers  is 
larger:  it  begins  in  the  transport  regime  and  extends  up  to  the  dispersive  region. 

The  aim  of  the  subsequent  MLE  simulations  is  to  check  the  GS  stability  predic¬ 
tions  obtained  in  Section  7.  To  this  end,  we  use  perturbed  GS  as  initial  conditions 
and  we  set  k  =  40  and  L  =  512  (2n)  to  approximate  an  infinite  geometry  by  means  of 
a  large  ring.  The  rest  of  the  parameters  are  left  unchanged,  ou2  —  1,  —  2,  and  the 

number  of  Fourier  modes  and  the  time  step  used  in  the  numerics  are  M-pourier  =  4096 
and  At  =  .01. 

The  first  three  plots  in  Fig.  38  corresponds  to  a  MLE  simulation  that  starts 
from  GS1  in  Fig.  26.  According  to  the  theoretical  results  from  Section  7,  this  is 
a  stable  soliton  and  this  is  also  what  is  obtained  from  the  MLE  simulation.  The 
time  evolution  of  ||£'||  (first  plot  in  Fig.  38)  exhibits  a  fast  oscillation  between  0 
(the  c  =  0  GS  have  equal  amplitudes)  and  its  constant  maximum  value.  And  the 
snapshots  of  the  spatial  profiles  of  E  at  two  different  times  (second  and  third  plot) 
show  that  the  GS  remains  practically  unchanged. 

The  initial  condition  used  in  the  last  two  plots  in  Fig.  38  is  the  GS2  in  Fig. 
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26.  This  GS  is  unstable  and  this  is  what  is  found  also  in  the  simulation.  The 
development  of  the  instability  produces  a  slow  time  modulation  in  the  fast  oscillation 
of  1 1 E 1 1  (fourth  plot)  and  the  distortion  of  the  GS  can  be  seen  in  the  snapshot  of  E 
(last  plot). 

The  GS2  and  GS4  in  Fig.  26  are  the  initial  conditions  of  the  simulations  presented 
in  Fig.  39.  The  GS2  is  a  stable  moving  soliton  and  this  is  confirmed  by  the  uniform 
oscillations  of  ||E'||  in  the  first  plot  of  Fig.  39  and  the  snapshots  of  E  (second  and 
third  plot)  that  show  a  translated  undistorted  GS.  Notice  that  the  GS  amplitudes 
are  not  equal  in  this  case  and  thus  the  minimum  of  fast  oscillation  of  ||i?||  is  not  0. 

The  GS4  is  unstable  and  again  this  is  reproduced  in  the  MLE  simulation:  a  slow 
time  oscillation  can  be  appreciated  in  the  time  evolution  of  ||i?||  (fourth  plot)  and 
the  MLE  solution  moves  away  from  the  GS  as  time  advances  (see  the  spatial  profile 
of  E  in  the  bottom  plot  of  Fig.  39). 

Finally,  in  agreement  with  the  stability  results  presented  in  Section  7,  the  MLE 
simulations  showed  no  sign  of  any  dispersive  instabilities  associated  with  GS.  The 
simulations  in  Figs.  38  and  39  correspond  to  negative  dispersion  {u  =  cu-),  but 
the  same  simulations  were  run  with  positive  dispersion  (c o  =  uj+)  and  no  significant 
qualitative  stability  differences  were  found. 
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Figure  33:  MLE  simulation  results  starting  from  a  CW  (k  =  1,  9  =  —  j  and  p2  =  1)  with  a  10-4 
perturbation.  From  top  to  bottom:  time  evolution  of  the  norm  of  E  for  u>~  and  u>+ ,  spatial  profiles 
of  E  at  t  =  75000  for  and  u+,  and  spatial  profile  of  E  at  t  =  160000  for  u>+  and  L  =  512-7T. 
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Figure  34:  MLE  simulation  results  starting  from  a  CW  (k  =  1,  9  =  —  j  and  p 2  =  4.5)  with  a  10  4 
perturbation.  Time  evolution  of  the  norm  of  E  for  u~  and  spatial  profile  of  E  at  t  =  88000. 


Figure  35:  MLE  simulation  results  starting  from  a  CW  (k  =  2,  9  =  —  \  and  p2  =  7.2)  with  a  10  4 
perturbation.  Time  evolution  of  the  norm  of  E  for  lo~  and  spatial  profile  of  E  at  t  =  4500. 
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Figure  36:  MLE  simulation  results  starting  from  a  CW  (k  =  1,  6  =  j  and  p2  =  1)  with  a  10-4 
perturbation.  From  top  to  bottom:  time  evolution  of  the  norm  of  E  for  and  to~ ,  spatial  profile 
of  E  at  t.  =  75000  and  u>+,  and  spatial  profile  of  E  at  t  =  30000  for  w-. 
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Figure  37:  MLE  simulation  results  starting  from  a  CW  (k  =  2,  9  =  j  and  p2  =  3)  with  a  10  4 
perturbation.  Time  evolution  of  the  norm  of  E  for  and  spatial  profile  of  E  at  t  =  18000. 
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Figure  38:  MLE  simulation  results  for  uj~  and  starting  from  a  GS  with  a  10-4  perturbation.  From 
top  to  bottom:  time  evolution  of  the  norm  of  E  starting  from  GS1  (c  =  0  and  7  =  cos  1),  spatial 
profiles  of  E  at  t  =  200  and  9000,  time  evolution  of  the  norm  of  E  starting  from  GS2  (c  =  0  and 
7  =  cos  2)  and  spatial  profile  of  E  at  t  =  2500. 


Figure  39:  MLE  simulation  results  for  u>  and  starting  from  a  GS  with  a  10  4  perturbation.  From 
top  to  bottom:  time  evolution  of  the  norm  of  E  starting  from  GS3  (c  =  0.5  and  7  =  -^cosl), 
spatial  profiles  of  E  at  t  =  150  and  6500,  time  evolution  of  the  norm  of  E  starting  from  GS4  (c  =  0 
and  7  =  cos  2)  and  spatial  profile  of  E  at  t  =  2000. 


9  Concluding  remarks 


Material  dispersion  terms  are  systematically  neglected  in  the  standard  amplitude 
equation  formulation  (NLCME)  used  in  the  literature  to  describe  the  weakly  non¬ 
linear  dynamics  of  light  propagation  in  fiber  gratings.  We  have  analyzed  in  this 
report  the  effect  of  dispersion  in  the  light  propagation  patterns  that  set  in  the  fiber 
grating.  To  be  more  precise,  this  work  focuses  on  the  effect  of  dispersion  on  the  sta¬ 
bility  of  two  important  families  of  solutions,  namely,  the  continuous  waves  (CW)  and 
the  gap  solitons  (GS).  We  have  obtained  precise  theoretical  linear  stability  results 
that  have  been  successfully  cross-checked  against  exhaustive  numerical  simulations 
of  the  amplitude  equations  with  small  dispersion  terms  (33)- (35)  and  of  the  com¬ 
plete  physical  model  equations  for  light  propagation  in  a  fiber  grating,  i.e.,  the  MLE 
(16)-(17).  The  development  of  the  numerical  codes  required  to  carry  out  the  two 
sets  of  simulations  described  above  has  also  been  part  of  the  present  research  effort. 

The  main  results  in  this  report  can  be  summarized  as  follows: 

•  The  stability  of  the  CW  is  drastically  affected  by  dispersion.  No  matter  how 
small  is  the  dispersion  or  what  sign  it  has,  there  are  always  stable  CW  according 
to  the  NLCME  formulation  that  are  dispersively  unstable.  The  destabilization 
produced  by  the  small  dispersive  terms  is  not  a  higher  order,  longer  time  effect; 
it  takes  place  in  the  time  scale  of  the  NLCME  and  the  associated  growth  rates 
remain  of  order  one  as  the  dispersion  coefficient  goes  to  zero. 

•  The  dispersion  induced  instability  generates  dispersive  scales  that  are  small 

(~  \/~L  1)  as  compared  with  typical  scale  of  the  NLCME  (~  L  1)  but 

still  large  as  compared  with  the  wavelength  of  the  basic  wavetrains  (~  1).  Once 
the  dispersive  scales  are  destabilized,  they  typically  spread  all  over  the  domain 
producing  a  very  complicated  spatio-temporal  dynamics  that  is  not  captured 
by  the  NLCME. 

•  This  behavior  is  the  result  of  the  competition  of  two  effects:  the  dominating 
advection  due  to  the  group  velocity  and  dispersion.  It  is  interesting  to  notice 
that  this  is  a  generic  situation  that  will  be  present  in  any  propagative  system 
with  spatial  reflection  symmetry  unless  we  manage  to  reduce  the  transport 
produced  by  the  group  velocity  (e.g.  by  tuning  the  parameters  to  approach  a 
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codimension  two  point). 

•  There  are  several  recent  rigorous  proofs  (see  [6]  and  [36])  that  establish  that 
the  solutions  of  the  NLCME  remain  asymptotically  close  to  solutions  of  the 
complete  physical  model  (the  MLE).  This  seems  to  be  in  contradiction  with 
the  results  of  this  report,  but  this  is  not  the  case.  These  proofs  do  not  contain 
any  stability  result  and  this  is  a  stability  issue:  when  a  stable  NLCME  solution 
is  dispersively  unstable,  its  close  MLE  solution  (without  dispersive  scales)  is 
also  unstable  and  the  MLE  dynamics  likes  to  move  away  from  this  solution 
and  develop  dispersive  scales,  as  the  amplitude  equations  with  small  dispersive 
terms  predict. 

•  The  GS  do  not  exhibit  any  dispersive  instability.  This  result  has  been  obtained 
from  a  linear  stability  analysis  and  has  been  checked  against  the  numerical  sim¬ 
ulations  of  the  amplitude  equations  with  dispersion  and  the  MLE.  The  reason 
for  this  behavior  lies  in  the  fact  that  this  is  a  propagative  instability  and  its 
growth  requires  the  coupling  that  comes  from  the  nonlinear  interaction  terms, 
but  these  terms  go  to  zero  as  the  perturbation  travels  away  from  the  soliton 
center  into  its  vanishing  tails. 

Summing  up  all  the  above  remarks,  we  can  conclude  that,  in  general,  the  standard 
NLCME  formulation  fails  to  predict  dynamics  of  the  system  and  the  material  dis¬ 
persion  terms  should  be  taken  into  account  (i.e.,  eqs.  (33)- (34)  should  be  used)  to 
appropriately  describe  the  weakly  nonlinear  regime  of  light  propagation  in  a  fiber 
grating. 
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Appendix  A.  Numerical  integration  of  the  ampli¬ 
tude  equations  (33)- (35). 

In  this  appendix  we  briefly  describe  the  numerical  method  used  to  integrate  the 
NLCME  with  dispersion  terms  (33)- (35). 

The  spatial  periodicity  of  the  problem  (33)- (35)  allow  us  to  expanded  its  solution 
in  finite  Fourier  series  [37,  38]  as 

k=-Of+ 1 

using  an  appropriately  large  number  of  Fourier  modes  N-p.  The  resulting  system  of 
ordinary  differential  equations  for  the  Fourier  modes  (ak  (t) ,  ak  (£))  can  be  written 
as 

da+ 

-jf  =  ctak  +  iKak  +  Ma5,  aJ ), 

^  =  ckak  +  iKat  +  MaJ ,  s+)’ 

where  =  ±i(2nk)  —  i£(2nk)2  and  the  nonlinear  terms  are  given  by 

nk(a+,dj)  =  [iA+(a\A+\2  +  |ff~|2)]fc  • 


In  order  to  avoid  numerical  instability  problems  associated  with  the  large  values 
of  the  linear  coefficients  for  large  k  [39],  we  first  rewrite  the  system  of  ODEs 
above  in  the  form 


d(e 

dt 

d(e-c7fq-) 

dt 


(inak  +  nk(a~j,  ))e  Cfct 
(inak  +nk(aj,a^))e~c'*t 


where  the  linear  diagonal  terms  are  integrated  exactly,  and  we  then  integrate  it  using 
a  standard  explicit  fourth  order  Runge-Kutta  method  [39]. The  resulting  integration 
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scheme  for  a  single  time  step  of  size  At  can  be  summarized  as  follows 


«*(*)  + 

At 

6 

[e*MK±  +  2e‘ 

^At/2K±  +  2ec 

Ktk  = 

/?( 

>■ fW), 

Kt  = 

/?( 

[ecfAt/2af{t)  + 

Kfk  = 

fkl 

ViA/  •-//  :(  /)  + 

— K ± ) 

2  2j)l 

K&  = 

/?( 

'ecfAtaf(t)  +  At  ecfAt'2K±), 

where  f^(a^)  =  ma +  nk(af,aj)  and  the  only  drawback  of  the  scheme  above  is 
the  fact  that  it  requires  to  store  the  arrays  of  coefficients  eck At/2,  which  take  into 
account  the  effect  of  the  linear  diagonal  terms. 

This  integration  procedure  always  handles  the  Fourier  representation  of  the  so¬ 
lution  The  solution  is  transformed  to  physical  space  only  during  the  com¬ 

putation  of  the  nonlinear  terms,  which  is  performed  in  the  following  steps 


l  ak  >  ak  . 


.  to  Phys. 


(A+,A~)^A±(a\A±\2  +  \A^\2) 


to  Fou. 


nk{aj  ,  a] 


where  the  products  of  the  amplitudes  are  computed  in  physical  space  and  the  so- 
called  2/3  rule  is  used  to  remove  the  aliasing  terms  (see  e.g.  [38]). 

This  numerical  procedure  has  been  implemented  in  a  C  code  and  the  FFTW 
subroutines  [40]  have  been  used  to  perform  the  Fourier  transforms.  The  typical 
resolutions  used  to  integrate  equations  (33)-(35)  with  dispersion  coefficient  in  the 
range  |e|  =  .001, . . . ,  .001/16  are  Np  =  256, . . . ,  1024  and  At  =  .01, . . . ,  .001. 

It  is  important  to  notice  that  the  number  of  Fourier  modes  used  for  the  numer¬ 
ical  integration  of  eqs.  (33)-(35)  should  not  be  too  large  since  otherwise  the  long 
wavelength  assumption  (19)  can  be  violated  [31].  The  correspondence  between  the 
wavenumber  of  the  Fourier  modes  of  the  solution  of  the  MLE  and  the  wavenumber 
of  the  modes  of  the  solution  of  the  amplitude  equations  (33)- (35)  is  given  by  (see 
expression  (13)) 

27T  1 

&MLE  =  ±1  +  -j-&AE,  with  (-)  ~  |£|  <C  1. 
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Figure  40:  Dispersion  relation  of  the  original  ML  equations  (solid  line)  and  of  the  amplitude 
equations  (33)-(35)  (dashed  line). 


And  the  parabolic  approximation  produced  by  the  amplitude  equations  to  the  true 
dispersion  relation  is  only  accurate  for  those  modes  of  the  MLE  with  wavenumbers 
near  &Mle  =  ±1  (see  the  dashed  line  in  Fig.  31).  In  order  to  avoid  spurious  solutions 
of  the  amplitude  equations  with  modes  corresponding  to  wavenumbers  away  from 
/cmle  =  ±1  (that  are  not  damped  because  of  the  absence  of  dissipation  in  the 
system),  and  to  ensure  that  the  dispersive  scales  (~  \f\s\)  are  properly  represented 
in  the  numerical  simulations,  the  following  relation  must  hold  between  the  dispersion 
coefficient  lei  -C  1  and  the  number  of  Fourier  modes  used  in  the  discretization  of  the 


amplitude  equations: 


The  above  problem  does  not  arise  in  strongly  dissipative  systems  in  which  the  dif¬ 
fusive  terms  wipe  out  the  high  wavenumber  modes. 
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Appendix  B.  Numerical 
(16)-(17). 

The  Maxwell-Lorentz  equations  (16)-(17) 
spatial  domain  of  size  L, 

E(x  +  L,t)  =  E(x ,  t ), 


integration  of  the  MLE 

with  periodic  boundary  conditions  in  a 

P{x  +  L,t)  =  P(x,t), 


are  integrated  using  in  a  numerical  procedure  completely  similar  to  the  one  described 
in  Appendix  A,  that  is,  using  spatial  discrete  Fourier  series  expansions  and  a  4th 
order  Runge-Kutta  temporal  scheme. 

In  order  to  do  this,  we  first  write  the  MLE  as  a  first  order  system  in  terms  of  E, 
E,  P  and  P, 


d  E 

~dt 

dE 

~dt 
d  P 

dt 

dP 

~dt 


Q2  t? 

-faj+uK  1  -  2Ancos(2x))P  -  u2p(nl  -  1  )E  -  WpP3, 

E, 

-cUp(l  -  2Ancos(2x))P  +  ujp(nl  -  1  )E  +  uj2P3, 

P, 


and  we  then  expand  the  solutions  in  discrete  Fourier  series, 


(E,  E,  P,  P)  =  et(t),Mt),Pk(t))ei(M,L)\ 

k 

using  an  appropriately  large  number  of  Fourier  modes.  Note  that  it  is  enough  to 
compute  the  modes  with  nonnegative  wavenumbers  because  the  electric  and  polar¬ 
ization  fields  are  real  fields  and  therefore  verify 


(^fcj  Pkj  Pk )  ip—ki  &—ki  P—ki  P—k )• 


The  Fourier  components  of  the  grating  term  in  the  system  above  can  be  explicitly 
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computed  to  give 


[2  cos(2a;)P]fc  =  [(ei2a;  +  e~i2x)  ^Pjeij(^x]k  = 

3 

l^2Pj(ei{j+nL)^)x  +  ei{j~nL)(^)x)}k  =  pk—nL  +pk+nL, 


where  n k  =  L/tt,  and  the  resulting  system  of  ODEs  for  the  Fourier  modes  of  the 
solution  can  be  finally  written  as 


dek 

dt 

dek 


=  ~(2irk/Lyek+uj2(pk  -  2A n(pk_nL  +  pk+nL))  ~  u;(n20  -  l)efc  -  c o2[P%, 


& k 5 


dt 

dpk  _ 

dpk 


=  -U2{pk  -  2A n(pk-nL  +pk+nL ))  +  0J2p{nl  -  l)ek  +  cu 2p[P% , 


dt 


=  Pk, 


with  k  —  0, . . . ,  M. 


This  system,  with  a  sufficiently  large  number  of  modes  M,  is  integrated  using  an 
explicit  4th  order  Runge-Kutta  method.  The  computation  of  the  nonlinear  term 
is  carried  out  in  physical  space  as  explained  in  Appendix  A  and  the  discrete  real 
Fourier  transforms  are  performed  using  the  FFTW  routines  [40]. 
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